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1.0  Justification 

Flutter  testing  is  both  hazardous  and  expensive,  but  it  is  required  to  verify  that  the  flutter 
boundary  lies  outside  of  the  operational  envelope.  It  is  imperative  that  this  testing  is  accomplished 
before  a  new  aircraft  design  is  released  to  an  unsuspecting  operational  pilot.  In  this  chapter,  the  history 
of  flutter  (and  other  aeroelastic  phenomenon)  will  be  reviewed  and  a  summary  of  the  popular  test 
techniques  will  be  presented.  Fiscal  factors  and  future  requirements  that  will  force  flutter  testing  to  be 
accomplished  more  efficiently  than  is  currently  possible  will  also  be  examined.  Finally  a  new  test 
methodology  that  utilizes  vibration  response  data  collected  during  a  longitudinal  acceleration  to  define 
the  frequency  and  damping  profile  within  the  operational  envelope  will  be  explored. 

The  objective  of  this  research  is  to  develop  a  test  methodology  that  will  reduce  test  cost  by 
simultaneously  reducing  the  number  of  flight  hours,  and  reducing  the  calendar  time  required  for 
accomplishing  flutter  testing. 

The  test  methodology  being  developed  is  designed  to  replace  the  traditional  flutter  clearance 
methodology  for  at  least  the  more  benign  parts  of  the  envelope.  Traditionally,  flutter  clearance  is 
accomplished  by  first  trimming  the  aircraft  at  a  predetermined  flight  condition.  A  forced  excitation  of 
the  aircraft  structure  is  then  initiated,  and  the  forced  inputs  and  responses  are  measured.  The  data  are 
analyzed  in  order  to  determine  the  stability  of  the  aircraft  structure,  and  a  decision  is  made  whether  or 
not  it  is  safe  to  continue  to  the  next  test  point.  This  process  is  expensive  and  time  consuming,  and  is 


the  motivation  for  this  research. 
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Following  is  a  short  introduction  of  flutter  and  traditional  flutter  test  techniques  that  is  presented 
in  order  to  better  draw  a  contrast  between  the  traditional  method  and  the  “accelerated  flutter  test” 
technique  that  is  the  subject  of  this  thesis. 


1.1  Description  of  Flutter 

Flutter  is  an  instability  caused  by  the  interaction  between  aerodynamic  forces  and  the  elastic 
structure. 

The  key  words/phrases  in  this  definition  are,  instability,  interaction,  aerodynamic  forces,  and 
elastic  structure. 

Starting  with  the  last  and  working  forward,  first  consider  the  elastic  structure.  The  linear  elastic 
structure  has  mass,  stiffness,  and  damping  so  you  can  think  of  the  system  as  an  oscillator.  The 
structure  is  often  represented  mathematically  as: 

Mx+Dx+Kx  =  0  (1.1) 

where  M  is  the  mass  matrix,  D  is  the  damping  matrix,  K  is  the  stiffness  matrix,  x  is  the  displacement 
vector,  and  a  dot  above  a  variable  represent  differentiation  with  respect  to  time.  This  is  also  known  as 
the  dry  structure,  or  the  structure  on  the  ground  (i.e.  without  airflow).  This  is  because  aerodynamic 
(“fluid”)  forces  are  not  acting  on  it. 

Next  consider  the  aerodynamic  forces.  That  is  the  force  that  the  air  exerts  on  the  structure  as  it 
rushes  past.  Mathematically  the  aerodynamic  force  can  be  considered  a  forcing  function  on  the  dry 
structure.  So  equation  (1.1)  becomes 
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M  x+Dx+Kx  =  F(t )  (1.2) 

Where  F(t)  is  the  time  dependent  aerodynamic  force,  and  t  is  the  variable  time. 

There  is  an  interaction  between  the  aerodynamic  forces  and  the  structure  because  the  aerodynamic 
force  is  dependent  on  the  motion  of  the  structure  and  the  motion  of  the  structure  is  driven  by  the 
aerodynamic  force.  Mathematically  the  aerodynamic  dependence  on  the  motion  of  the  structure  can  be 
expressed  as 


F(t)  =  XFM  x+XFd  x+X FKx  +  Fr(t)  (1.3) 

Where  X  is  the  dynamic  pressure,  FM  ,  Fd  , and  Fk  ,  are  matrices  representing  aerodynamic 
operators,  and  Fr  represents  unsteady  aerodynamic  forces  which  are  independent  of  the  aircraft 
motion.  An  example  of  a  force  that  would  be  included  in  FT  is  the  unsteady  force  due  to  pressure 
fluctuations  induced  by  high  angle  of  attack.  Combining  equations  (1.2)  and  (1.3)  yields: 


M  x+  D  x+  Kx  =  XFm  XFd  x+  XFkx  +  FT  (t)  (1.4) 

Next  move  the  terms  that  are  dependent  on  the  motion  of  the  structure  to  the  left  hand  side  of  the 
equation. 


(M-XFm)x+(D-AFd)x+(K-~XFk)x  =  FT(t)  (1.5) 

In  equation  (1.5)  it  is  easy  to  see  the  dynamic  pressure  (speed  of  the  aircraft)  plays  a  significant 
role  in  the  mass,  damping,  and  stiffness  of  the  system.  This  is  what  is  called  the  wet  structure. 
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Furthermore  it  is  easy  to  imagine  a  case  where  damping  for  at  least  one  mode  in  this  system  could 
be  driven  to  instability  by  the  aerodynamic  force  acting  on  the  structure.  In  the  case  where  this 
instability  manifests  itself  the  aircraft  is  said  to  flutter.  For  more  information  on  classical  flutter 
development  see  references  [2,28,29,31]. 

In  order  to  illuminate  the  differences  between  traditional  flutter  testing  and  the  new  method  being 
examined  in  this  work,  this  chapter  will  present  a  brief  history  of  flutter  (in  section  1 .2),  review  the 
traditional  flutter  testing  methodology  and  cost  (in  sections  1.3,  and  1.4  respectively),  and  explore 
manufacturing  developments  that  will  make  it  difficult  to  utilize  the  traditional  approach  in  the  future 
(in  section  1.5).  Finally  this  chapter  will  present  the  objectives  of  this  research  (in  section  1.6). 


1.2  A  History  of  Flutter 

Aeroelastic  problems  date  back  as  far  as  four  centuries  ago,  to  the  windmills  in  Holland.  The 
spars  that  were  located  at  the  mid-chord  were  moved  to  the  quarter-chord  to  prevent  twisting.  Also 
some  19th  century  bridges  were  weak  in  torsion  and  were  destroyed  by  aerodynamic  effects.  For  more 
on  the  history  of  flutter  see  [1]. 

The  first  occurrence  of  an  aeroelastic  instability  of  an  aircraft  was  actually  before  the  Wright 
brothers  first  flew.  Samuel  Langley  launched  his  monoplane  from  a  houseboat  in  the  Potomac  River 
and  suffered  a  structural  failure  due  to  divergence.  Although  divergence  is  not  flutter,  it  is  closely 
related.  In  fact  a  flutter  analysis  will  give  the  divergence  speed  at  zero  frequency.  So  divergence  can 
be  thought  of  as  zero  frequency  flutter. 


During  World  War  I  the  Handley  Page  bomber  had  a  control  surface  flutter  problem  that  was 
described  as  a  coupling  between  the  aft  fuselage  and  antisymetric  elevator  motion.  As  a  result  of  this 


and  similar  problems  of  the  DH-9  in  1917  aircraft  designs  began  to  incorporate  a  stiff  torque  tube 
connecting  the  elevators  so  the  frequency  of  the  antisymmetric  elevator  mode  would  be  much  higher 
than  the  torsional  mode  for  the  fuselage. 

In  1925  two  aircraft  designs,  the  Gloster  Grebe  and  the  Gloster  Gamecock,  were  found  to  be 
unstable  due  to  wing  aileron  flutter.  The  problem  was  fixed  by  changing  the  shape  of  the  aileron. 
Effectively  the  mass  of  the  aileron  was  moved  forward.  More  recently  using  balance  weights  to  move 
the  center  of  gravity  of  the  aileron  has  alleviated  this  type  of  flutter  simply. 

The  first  use  of  a  scaled  model  to  determine  the  critical  flutter  speed  (the  lowest  speed  at  which 
flutter  onset  can  occur)  was  1928.  A  one  third  scale  model  with  the  same  mass  density  and  one  ninth 
the  stiffness  was  used,  and  the  resulting  flutter  speeds  and  frequencies  of  the  model  correlated  well 
with  those  of  the  full  scale  aircraft. 

Many  flutter  incidences  in  the  late  1920s  and  early  1930s  involved  air  racers.  One  account,  in 
1934,  the  pilot  repeatedly  encountered  wing-tip  flutter.  Each  time  he  reduced  the  wingspan  by  cutting 
off  the  tip  of  the  wing.  Eventually  the  wing  area  was  reduced  from  78  square  feet  to  42  square  feet, 
but  the  airplane  was  then  flutter-free. 

Von  Schlippe  was  the  first  to  accomplish  systematic  flutter  flight-testing  in  1935.  Prior  to  this 
time  flutter  flight-testing  consisted  of  diving  the  aircraft  to  its  maximum  velocity  and  hoping  for  the 
best.  By  the  late  1940s  systematic  flight  flutter  testing  had  gained  wide  acceptance  in  the  industry. 

The  first  attempt  to  use  a  flutter  model  to  simulate  free  flight  conditions  was  attempted  in  1944. 

The  first  example  of  panel  flutter  occurred  during  development  of  the  V-2  rocket  between  1942 
and  1944.  Panel  flutter  depends  on  many  parameters,  which  include  Mach  number  and  boundary 
layer,  but  is  especially  sensitive  to  any  thermal  or  compressive  effects  that  tend  to  create  local  buckling 


in  the  skin. 
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During  flight-testing  of  the  P-80  aircraft,  in  1944,  a  type  of  flutter  called  aileron  buzz  was 
encountered.  This  phenomenon  is  a  single-degree-of-ffeedom  form  of  flutter  caused  by  the  coupling 
of  aileron  rotation  and  chordwise  motion  of  shock  waves  on  the  wing. 

Propeller  whirl  flutter  is  suspected  in  two  fatal  crashes  of  the  Lockheed  Electra  turboprop 
transports.  Propeller  whirl  flutter  involves  the  gyroscopic  precession  of  a  flexibly  mounted  engine- 
propeller  system.  It  was  first  explored  academically  in  1938,  but  was  thought  to  be  only  of  academic 
interest  until  the  Electra  crashes  in  1959. 


1.3  Flutter  Test  Techniques 

Flutter  testing  can  be  decomposed  into  four  areas;  1)  excitation  type,  2)  excitation  method,  3) 
response  measurement,  and  4)  data  analysis.  The  capabilities  and  limitations  of  the  techniques  used  to 
accomplish  these  four  areas  have  driven  flutter  flight  test  engineers  to  the  usual  flutter  testing 
techniques  that  have  been  used  for  decades. 

The  usual  method  used  for  flutter  testing  aircraft  is  to  fly  stabilized  test  points  at  predetermined 
conditions.  For  example,  to  clear  the  aircraft  at  10,000  feet  for  a  given  configuration,  it  may  be 
necessary  to  conduct  stabilized  test  points  at  eight  to  twelve  speeds  at  that  altitude.  After  the  first  two 
or  three  test  points  are  accomplished  within  the  analytically  cleared  envelope,  the  damping  trend  is 
extrapolated  and  a  determination  is  made  whether  or  not  it  is  safe  to  continue  to  the  next  test  point.  If 
it  is  safe  to  go  to  the  next  test  point,  the  process  is  continued  until  the  top  speed  for  that  altitude  is 
achieved.  It  is  easy  to  see,  if  several  configurations  have  to  be  cleared  at  several  altitudes  the  full  test 
matrix  becomes  very  time  consuming  and  expensive.  For  more  information  on  flutter  flight  test 
procedures  see  references  [24,26,30,36]. 
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1.3.1  Excitation  Type 

Excitation  type  describes  the  input  signal  that  generates  the  vibration  response  of  the  airframe. 

The  excitation  types  that  are  commonly  used  are  decay,  impulse,  harmonic  excitation  (sine  dwell),  sine 
sweep,  tailored  inputs,  and  random  inputs.  For  more  information  on  excitation  types  see  references 
[37,38,39]. 

1.3. 1.1  Decay  Data 

Decay  data  are  what  you  get  in  the  absence  of  an  excitation  signal.  Typically  the  aircraft  structure 
is  excited  by  some  other  type  of  excitation  after  which  the  excitation  is  abruptly  halted.  The  data  are 
recorded  and  analyzed  for  the  period  after  the  excitation  has  been  halted. 

1 .3. 1 .2  Impulse  Excitations 

Impulse  excitations  are  typically  used  as  a  method  to  generate  decay  data.  Striking  the  control 
stick  with  a  hand  or  a  mallet,  kicking  the  rudder,  or  detonating  pyrotechnic  charges  located  on  the  wing 
of  the  aircraft,  can  produce  the  impulse.  Striking  the  control  stick  with  a  hand  or  kicking  the  rudder 
cannot  adequately  excite  structural  modes  with  frequencies  higher  than  about  5-6  Hz.  Figure  1 . 1 
shows  a  typical  impulse  excitation  along  with  the  decay  response. 
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Figure  1.1:  Typical  impulse  excitation  and  decay  response 


1 .3 . 1.3  Harmonic  Excitations 

Harmonic  excitation  (sine  dwell)  is  a  sustained  sinusoidal  signal  that  is  tuned  to  the  harmonic 
frequency  to  be  excited.  This  excitation  is  useful  because  it  typically  excites  only  one  mode,  and  it 
puts  more  energy  into  that  one  mode.  Unlike  decay  data  and  impulse  data,  harmonic  data  are  typically 
recorded  during  the  excitation  and  therefore  the  analysis  methods  should  account  for  the  driving  force. 

1.3. 1.4  Sine  Sweep 

Sine  sweep  is  similar  to  sine  dwell  except  the  frequency  of  the  sine  wave  is  varied  with  time. 
Sometimes  the  variation  is  linear  and  sometimes  the  variation  is  exponential.  The  exponential  sine 
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sweep  is  normally  preferred  because  it  puts  an  equal  amount  of  energy  into  each  frequency.  Figure  1.2 
shows  a  typical  sine  sweep.  In  this  case  the  frequency  varies  linearly  with  time. 


Figure  1.2:  Typical  sine  sweep  using  linear  variation  in  frequency. 


1 .3. 1 .5  Tailored  Inputs 

Tailored  inputs  focus  energy  on  predetermined  frequency  bands.  These  are  of  academic  interest, 
however,  they  are  not  used  much  in  practice  because  they  require  knowing  the  frequencies  of  the 
system  under  test  in  advance. 

1 .3 . 1 .6  Random  Inputs 

Random  inputs  use  a  pseudo-random  signal  to  excite  the  structure  of  the  aircraft.  This  can  be  a 
very  quick  way  to  identify  the  structure,  however  the  energy  is  spread  across  a  wide  spectrum  and  may 
not  adequately  excite  modes. 
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1 .3.2  Excitation  Methods 

Excitation  methods  that  are  commonly  used  fall  into  three  categories;  1)  aerodynamic  methods,  2) 
moving  mass  methods,  and  3)  pyrotechnic  methods.  The  aerodynamic  methods  are  control  surface 
excitation,  oscillating  vane  excitation,  and  random  turbulence  excitation.  The  moving  mass  methods 
include  inertial  excitation  and  electrodynamic  excitation.  Pyrotechnics  methods  use  small  explosive 
charges  to  generate  an  excitation.  Sometimes  these  pyrotechnic  devices  are  called  bonkers.  For  more 
information  on  excitation  methods  see  references  [37,38,39]. 

1.3. 2.1  Aerodynamic  Excitation  Methods 

The  first  of  the  aerodynamic  methods  used  to  excite  the  aircraft  structure  is  control  surface 
excitations.  Control  surface  excitations  can  be  pilot  induced  or  flutter  excitation  system  induced.  Pilot 
induced  excitations  are  low  quality,  unrepeatable,  and  normally  contain  only  frequencies  below  5-6 
Hz.  Flutter  excitation  systems  are  capable  of  producing  repeatable  signals  with  any  frequency  content 
the  aircraft  is  capable  of  producing. 

Oscillating  vane  excitations  can  be  produced  by  an  airfoil  or  slotted  tubes.  The  airfoil  can 
produce  any  type  of  excitation,  but  the  slotted  tubes  can  only  produce  sine  dwell  and  sine  sweep 
excitations.  The  slotted  tubes  also  have  a  tendency  to  get  stuck  due  to  aerodynamic  stresses  in  the 
tubes.  Figure  1.3  shows  a  typical  aerodynamic  vane  excitation  system  that  utilizes  airfoils  to  produce 


an  excitation  force. 
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Figure  1.3:  Aerodynamic  excitation  system  utilizing  oscillating 
airfoils  to  produce  an  excitation  force 


Random  turbulence  excitation  is  the  polite  name  for  flying  around  looking  for  gusts.  This  method 
is  only  used  when  there  is  no  alternative  (i.e.  the  program  cannot  afford  to  properly  instrument  the 
aircraft,  or  when  the  aircraft  excitation  system  cannot  excite  the  frequency  of  concern).  The  signal  to 
noise  ratio  is  very  low  for  this  type  of  data  so  the  data  have  to  be  averaged.  Care  must  be  taken  to 
assure  that  the  data  that  are  averaged  is  in  phase.  This  can  be  done  using  an  algorithm  called  random 
decrement  if  there  are  only  one  or  two  dominate  modes  in  the  data.  The  random  decrement  algorithm 
is  basically  a  triggering  technique.  Data  are  collected  from  the  time  the  trigger  is  tripped  for  a 
predetermined  number  of  samples.  Filters  are  sometimes  used  to  assure  there  will  only  be  one  or  two 
modes  in  the  data.  After  the  resulting  data  are  averaged  it  is  treated  like  decay  data. 

1 . 3 .2 .2  Moving  Mass  Excitation  Methods 


Moving  mass  methods  include  inertial  excitation,  another  name  for  an  unbalanced  mass.  An 
unbalanced  mass  spinning  on  a  shaft  creates  an  oscillating  force  on  the  structure.  Unless  the  mass  can 
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be  moved  the  force  increases  as  the  shaft  rotates  faster.  This  results  in  a  system  that  has  a  small 
frequency  range  in  which  it  is  effective.  This  is  typically  a  low  budget  approach  to  flight-testing.  If 
the  program  can  afford  it,  it  is  better  to  use  one  of  the  aerodynamic  methods  discussed  above  with  the 
exception  of  random  turbulence. 

Electrodynamic  excitation  is  similar  to  the  unbalanced  mass  except  the  mass  is  driven  by  an 
electrodynamic  shaker.  This  is  the  excitation  method  that  was  used  on  the  B-1B  flutter  envelope 
expansion.  Figure  1 .4  shows  a  cartoon  of  the  inertial  shaker  system  that  was  used  on  the  B-1B.  The 
benefits  of  using  this  method  include  good  frequency  range,  measurable  input  force,  and  repeatable 
execution.  The  primary  drawback  is  the  mass  required  to  generate  a  suitable  force  on  the  structure  can 
change  the  modal  characteristics. 


hydmulicaty  acluafad  movaW#  mtn  (wand) 


Figure  1.4:  Inertial  excitation  system  used  on  the  B-1B 
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1 .3.2.3  Pyrotechnic  Excitation  Methods 

Bonkers  are  not  used  in  the  United  States  because  of  Safety  of  Flight  considerations.  However, 
Europe  has  a  slightly  different  philosophy  that  allows  the  use  of  pyrotechnics.  The  bonkers  produce  a 
nice  impulse  force  with  good  frequency  content,  and  they  can  be  detonated  together  to  excite  a 
particular  mode.  The  only  problems  with  the  bonkers  are  the  safety  implications  and  the  number  of 
test  points  that  can  be  accomplished  in  one  mission  is  determined  by  the  number  of  bonkers  that  can  be 
attached  to  the  wings. 

1.3.3  Response  Measurement 

There  are  only  two  kinds  of  sensors  that  are  typically  used  to  measure  the  response  of  the  system. 
Accelerometers  are  used  wherever  the  displacements  tend  to  be  large  (e.g.  wing  tips),  and  strain  gauges 
are  used  where  the  strain  tends  to  be  large  (e.g.  wing  root).  Experience  suggests  the  accelerometer  data 
are  the  best  suited  for  identifying  frequency  and  damping.  Typical  sample  rates  for  flutter  data  range 
from  100  Hertz  to  800  Hertz.  For  more  information  on  response  measurement  and  how  the  data  is 
transmitted  to  the  control  room  see  references  [32,33,34,35]. 

1.3.4  Data  Analysis 

There  are  several  techniques  widely  used  for  analyzing  flutter  test  data.  The  selection  of  analysis 
technique  is  highly  dependent  on  the  excitation  method  that  is  used.  Some  of  the  common  data 
analysis  techniques  are  half  power,  time  history  curve  fit,  frequency  curve  fit,  random  decrement, 
pseudo-random  decrement,  and  log  decrement.  For  more  information  on  data  analysis  methods  see 
references  [40,41,42]. 


1.3. 4.1  Half  Power  Analysis 
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Half  power  uses  the  width  of  the  modal  lobe  at  the  half  power  line  of  the  Power  Spectral  Density 
(PSD),  the  auto-spectrum  or  the  Fast  Fourier  Transform  (FFT),  to  determine  the  damping  value  for  that 
particular  mode.  Figure  1.5  shows  how  the  half  power  line  is  defined  for  the  power  spectral  density  or 
an  auto-spectrum  (cases  where  the  amplitude  has  squared  units).  Figure  1.6  shows  how  the  half  power 
line  is  defined  for  an  FFT  (cases  where  the  amplitude  does  not  have  squared  units).  Once  the  lobe 
width  at  the  half  power  line  is  determined  as  in  equation  ()  then  the  damping  can  be  determined  using 
equation  ().  This  technique  is  highly  dependent  on  the  analyst  and  it  becomes  significantly  less 
reliable  when  considering  closely  spaced  modes.  This  is  particularly  problematic  because  two  closely 
spaced  modes  are  often  an  indication  of  impending  flutter  conditions.  The  advantage  of  the  half  power 
analysis  technique  is,  it  can  be  used  with  most  data  collection  methods.  However  this  flexibility  comes 
at  the  cost  of  confidence  in  the  results. 


Figure  1.5:  Half  power  method  for  Power  Spectral  Densities  and  Auto 
spectrum. 


Magnitude 
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Frequency 

Figure  1.6:  Half  power  method  for  Fast  Fourier  Transforms. 


&F  =  Fhi-Fh=co 2~a\  (1-6) 


C  = 


2  2 

co2  -cox 
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(1.7) 


1 .3 .4.2  Log  Decrement  Analysis 

Log  decrement  uses  decay  data  to  assess  the  damping  of  the  mode  of  interest  taking  advantage  of 
the  relationship  between  the  damping  value  and  the  number  of  cycles  ( N )  it  takes  for  the  amplitude  of 
the  mode  at  time  zero  (  Aq  )  to  be  reduced  to  some  fractional  value  ( AN ).  Equation  (1.8)  gives  the 
algorithm  that  is  used  to  determine  damping  from  the  decay  data.  Sometimes,  in  order  to  expedite  the 
estimation  of  damping,  equation  (1.9)  is  used.  In  equation  (1.9)  N1/2  is  the  number  of  cycles  it  takes 
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to  reduce  the  amplitude  to  Aq /2 .  Sometimes  band  pass  filtering  is  used  to  isolate  modes.  Log 

decrement  offers  consistent  results,  and  is  much  less  dependent  on  the  analyst  than  half  power. 
However,  like  half  power,  it  is  not  very  effective  when  dealing  with  closely  spaced  modes.  This 
method  is  also  limited  to  decay  data,  which  does  not  always  contain  the  modal  energy  necessary  to 
achieve  accurate  results. 


C  =  HA/An)/(2ttN)  (1.8) 


<Z  =  0.n/NV2  (1.9) 

1 .3.4.3  Time  History  Curve  Fit  Analysis 

Time  history  curve  fit  is  a  least  squares  curve  fit  of  the  time  history  data  that  assumes  an 
underlying  function  comprised  of  damped  sinusoids  (Equation  (1.10)  gives  the  general  form  of  the 
assumed  function  for  time  history  curve  fit). 

y(t)  =  c + X sin(fJ +(P„)  (iio) 

n=l 


Where: 

4,  = 

maximum  amplitude  of  mode  n 

c 

zero  offset  correction 

t 

time  in  seconds 

y(t)  = 

value  of  signal  at  time  t 

N 

maximum  mode  number 

<Pn 

phase  angle  of  mode  n  in  radians 
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natural  frequency  of  mode  n  in  radians/second,  and 
damping  ratio  of  mode  n 


The  method  iteratively  improves  the  damping  estimate  and  then  the  frequency  estimate.  The 
method  requires  an  a  priory  estimate  of  the  frequency  which  is  normally  accomplished  by  using  a  Fast 
Fourier  Transform  of  the  flight  test  data.  Because  it  assumes  damped  sinusoids  this  method  is  only 
useful  for  decay  data.  When  this  method  converges  to  a  solution  it  provides  consistent  and  accurate 
results. 


1 .3.4.4  Frequency  Curve  Fit  Analysis 


Frequency  curve  fit  is  a  least  squares  curve  fit  of  the  cross  spectrum  of  the  flight  test  data  that 


assumes  an  underlying  function  that  is  the  ratio  of  complex  polynomials  which  is  shown  in  equation  (). 


G, .  v  _  Po+PiUa>)  +  P2U*>)2+  —  =  P(ja) 

\  +  qx(jco)  +  q2{jcof  +■■■  Q(jco) 

C Po  -  Pi  +  P± (ft>)4  •  •  •)  +  ja(Pi  -  ft  O)2  +  Ps  O)4  •  •  •) 

(1  -  q2  (a)2  +  q4  (co)A  •  •  •)  +  jco(qx  -  q3  (a)2  +  q5  (co)4  •  •  •) 

=  a  +  JC0/]  *  F(joj)  =  R(ca)  +  jI(G>) 

(7  +  JCOT 


Where: 


G{jco) 

p's 


^approximating  function  representing  the  transfer  function 
^unknown  numerator  coefficients  to  be  solved  for 


f 
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q '  S  =unknown  denominator  coefficients  to  be  solved  for 

F(jco)  =cross  spectrum  of  the  flight  test  data  which  approximates  the  transfer  function 

R(a>)  =real  part  of  F(jo)) 

1(G))  imaginary  part  of  F(jG)) 

The  damping  and  frequency  values  can  be  calculated  from  the  roots  of  the  denominator,  which 
give  the  poles  of  the  system.  This  method  requires  the  interaction  of  the  analyst,  and  two  to  four 
minutes  worth  of  data.  This  method  is  only  useful  when  an  excitation  system  is  used  to  excite  the 
structure.  For  more  on  this  analysis  technique  see  references  [22,23]. 

1.3. 4.5  Random  Decrement  Analysis 

Random  decrement  is  not  an  analysis  technique  as  much  as  it  is  a  preconditioner.  This  technique 
is  used  on  random  turbulence  data  in  order  to  average  out  the  unwanted  signal.  It  works  by  taking  a 
block  of  data  starting  at  the  point  the  amplitude  of  the  data  crosses  a  trigger  value  that  is  set  by  the 
analyst.  In  the  end  all  of  the  blocks  of  data  are  averaged.  The  resulting  signal  is  referred  to  as  the 
randomdec  signature.  The  randomdec  signature  can  then  be  used  with  any  analysis  technique  that 
utilizes  decay  data.  Typically  this  technique  is  used  in  conjunction  with  the  log  decrement  method. 
This  technique  typically  requires  several  minutes  of  flight  test  data  to  achieve  results.  For  more 
detailed  information  regarding  random  decrement  analysis  see  reference  [25]. 

1 .3.4.6  Psudo-Random  Decrement  Analysis 

Psudo-random  decrement  is  much  like  the  random  decrement  method.  In  this  case  an  excitation 
system  is  used  to  induce  vibration  into  the  structure,  then  the  excitation  system  is  abruptly  halted. 
Instead  of  using  a  trigger  level  to  determine  when  the  blocks  of  data  are  started,  this  technique  starts 
the  block  of  data  a  few  samples  after  the  excitation  system  is  halted.  This  method  can  identify  a  single 
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mode  of  interest  using  data  from  a  few  excitations.  The  total  elapsed  time  necessary  to  collect  the  data 
from  these  excitations  is  about  15  to  30  seconds.  This  is  a  very  useful  method  as  long  as  there  are  only 
one  or  two  modes  of  interest. 


1.4  Cost  of  Flutter  Testing 

The  F-16  Block  60  program  will  be  used  as  an  example  of  the  cost  of  flight  testing,  and  the 
potential  benefit  of  using  the  flight  test  and  simulation  techniques  described  in  this  thesis.  The  F-16 
Block  60  is  a  significant  modification  to  the  current  F-16  that  is  being  produced  by  Lockheed  Aircraft 
for  sale  to  the  United  Arab  Emirates.  Structurally  the  Block  60  will  not  change  much,  however,  the 
surface  geometry  will  change  somewhat  in  order  to  increase  the  fuel  capacity.  Barring  unexpected 
complications,  flutter  clearance  will  be  straightforward.  However,  this  program  intends  to  certify  1500 
weapons  configurations  including  fully  loaded  configurations,  as  well  as  configurations  that  occur  as  a 
result  of  selectively  dispensing  weapons.  Using  current  methods,  only  the  10%,  or  150  loads  that  are 
both  stable  and  easily  simulated,  can  be  cleared  using  simulation.  For  example  low  subsonic 
conditions  with  lightweight  and  stable  stores.  The  rest  (1350  loads)  should  be  cleared  through  flight 
test.  I  say  should  be  cleared  through  flight  test  because  it  is  prohibitive  to  do  this  much  flutter 
clearance  work  for  a  single  program.  In  reality  what  would  happen  is  the  program  will  assume  more 
risk  by  choosing  what  is  believed  to  be  the  worst  case  loads  and  leaving  many  loads  untested  even 
though  these  loads  are  potentially  hazardous.  Therefore  the  benefits  improving  the  efficiency  of  the 
flutter  clearance  process  will  be  a  combination  of  cost  avoidance  and  risk  reduction. 

The  analysis  presented  below  is  an  attempt  to  quantify  the  benefits. 

Starting  with  the  requirements  to  test  without  clearing  the  configuration  using  accelerating  flight- 
test  data,  or  coupled  field  analysis,  1350  configurations  must  be  cleared  using  flight  test.  Normally  it 
would  take  at  least  3  hours  to  clear  one  configuration.  This  analysis  will  use  the  condition  that  each 
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configuration  is  cleared,  on  average,  in  2  hours.  Therefore  testing  will  take  about  2700  hours  of  flight 
time. 


Traditionally,  predictive  flutter  clearance  analysis  is  done  using  a  structural  dynamics,  and 
aerodynamics  models  that  are  very  low  fidelity.  It  has  been  necessary  in  the  past  to  use  these  low 
fidelity  models  due  to  computational  power  that  was  available  at  the  time.  Supposing  that  the  cost  of 
doing  this  traditional  analysis  is  negligible.  It  is  safe  to  say  that  the  total  cost  of  the  flight  test  program 
done  to  completion  would  be  at  least  $15M. 

The  lowest  rate  the  Air  Force  Flight  Test  Center  would  charge  for  flight  testing  an  F-16  aircraft  is 
about  $5600  per  hour  (as  of  2001).  This  rate  does  not  include  the  control  room,  range  time,  telemetry, 
or  a  chase  aircraft.  All  of  which  are  required  for  hazardous  flight-testing.  The  cost  of  2700  hours  at 
$5600  per  hour  is  just  over  $15M.  What  is  worse  than  the  cost,  is  the  schedule  impact  of  doing  that 
much  testing.  If  eight  flights  per  week  are  flown,  and  one  hour  of  testing  is  accomplished  during  each 
flight  it  will  take  330  weeks  (over  6  years)  to  complete  the  flight  test  program. 

Compare  the  program  cost  and  schedule  above  to  that  of  using  coupled  field  analysis  and  data 
from  accelerated  flight.  If  the  fidelity  of  coupled  field  analysis  is  such  that  Limit  Cycle  Oscillation  can 
be  dependably  predicted  the  total  number  of  flight  test  points  that  must  be  flight-tested  could  be  cut  by 
an  order  of  magnitude  (or  90%).  That  means  that  150  configuration,  would  have  to  be  flight-tested. 
Also,  using  data  from  accelerated  flight  could  reduce  the  number  of  hours  to  clear  a  configuration  to 
one.  Therefore  it  would  take  150  flight  hours  to  complete  the  flight-testing. 

Realizing  that  the  complexity  per  flight  hour  has  increased,  allowance  needs  to  be  made  for  the 
rate  to  increase.  So,  for  this  case  the  rate  of  $  1 0,000/hour  will  be  used.  Therefore  the  total  cost  of  the 
flight  test  would  be  $  1 .5M.  Add  to  that  $  1M  for  conducting  the  simulations.  The  total  cost  of  the 
program  is  $2.5M.  The  total  return  would  be  at  least  $  12.5M. 
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As  far  as  schedule  goes,  a  more  realistic  sortie  rate  of  three  flights  per  week  will  still  complete  the 
program  in  50  weeks,  or  less  than  one  year.  With  these  cost  and  schedule  savings,  the  program  could 
truly  clear  the  entire  1500  configurations,  and  not  risk  the  lives  of  operational  users,  and  the  program 
could  still  come  in  under  cost  and  schedule. 


1.5  New  Technologies  and  Tactics 

New  technologies  such  as  flexible  manufacturing,  rapid  prototyping,  and  powder  metallurgy  will 
make  it  practical  to  build  aircraft  without  going  into  large-scale  production.  For  example  if  you  want 
to  build  one  or  two  of  a  kind  aircraft  for  a  specific  mission  it  will  be  possible  to  do  that.  This  is  also 
consistent  with  some  of  the  conceptual  tactics  that  are  starting  to  surface.  It  is  feasible  to  counter  a 
handful  of  fighter  designs  that  are  currently  in  the  US  arsenal,  but  if  hundreds  of  fighter  designs  are 
available,  counter  measures  become  much  less  tractable.  All  of  this  pretty  much  assures  there  will  be 
an  increased  need  to  do  flutter  testing  in  the  future.  With  a  shrinking  flight  test  force  it  will  become 
critical  that  clearance  of  the  flutter  envelope  for  new  aircraft  be  accomplished  more  efficiently  than 
ever  before.  Flutter  test  techniques  that  are  currently  used  differ  only  slightly  from  the  techniques  that 
were  used  50  years  ago.  Computer  technology  has  been  used  to  expedite  the  data  analysis  process,  but 
has  not  been  used  to  expedite  the  flight  test  process.  This  research  is  focused  on  using  high  fidelity 
simulations  in  conjunction  with  a  new  flight  test  technique  to  demonstrate  a  revolutionary  flight  test 
approach  capable  of  achieving  the  increased  efficiency  that  the  future  demands. 


1.6  Objectives 

In  the  new  test  methodology  the  aircraft  would  stabilize  at  a  low  airspeed  and  predetermined 
altitude/configuration.  While  the  aircraft  maintains  altitude,  an  excitation  system  would  be  turned  on 
and  the  aircraft  would  begin  to  accelerate.  As  the  airspeed  increases  the  frequencies  and  damping 
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values  would  continually  be  calculated  and  displayed  very  near  real  time.  The  aircraft  would  continue 
to  accelerate  until  it  reaches  a  predetermined  velocity  as  long  as  the  analysis  indicates  it  is  safe  to 
continue.  There  are  two  big  advantages  to  testing  using  this  methodology  rather  than  traditional  flutter 
test  methods.  The  first  is  simply  test  efficiency.  One  altitude/configuration  combination  could 
potentially  be  cleared  in  a  matter  of  minutes  rather  than  the  hours  it  can  take  in  some  cases.  The 
second  is  the  resolution  of  data.  In  the  traditional  approach,  damping  estimates  are  only  available  for 
the  test  point  conditions  that  have  been  accomplished.  In  the  accelerated  approach  there  is  a  potential 
that  the  frequency  resolution  of  the  damping  curve  would  be  almost  continuous. 

In  order  to  accomplish  this  objective  there  are  several  issues  that  must  be  addressed. 

1)  The  first  issue  is  whether  or  not  the  frequencies  and  dampings  are  affected  by 
the  acceleration  of  the  aircraft. 

2)  The  second  issue  is  whether  or  not  parameter  identification  algorithms  can 
extract  frequency  and  damping  values  from  the  time  varying  data. 

3)  The  third  issue  is  the  need  to  anticipate  the  vibration  response  at  airspeeds  (or 
Mach  numbers)  beyond  which  the  aircraft  has  accelerated. 

4)  The  fourth  issue  is  developing  formal  criteria  to  determine  when  the  aircraft 
needs  to  end  the  acceleration  and  terminate  the  test  point. 

Before  these  issues  are  examined  in  detail,  there  are  certain  requirements  that  must  be  observed  if 
this  research  is  to  become  a  usable  flight  test  tool.  These  requirements  will  be  examined  in  chapter  2. 

The  first  of  these  issues  will  be  addressed  in  chapter  3.  The  typical  wing  section  will  be 
introduced  in  this  chapter,  and  the  effects  of  acceleration  on  the  typical  wing  section  will  be  explored 
analytically.  Using  the  typical  wing  section,  a  model  will  be  developed  that  will  approximate  the  F- 16 
fighter  aircraft.  Finally  it  will  be  shown  that,  for  this  practical  case,  the  effects  of  the  accelerations  are 
negligible. 
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The  second  of  these  issues  will  be  addressed  in  chapter  4.  This  chapter  will  introduce  the  idea  of 
windowing  and  its  effect  on  time  varying  data.  A  summary  will  be  give  of  the  Eiganvector  Realization 
Algorithm  (ERA),  and  its  variant  FastERA  and  the  utility  of  FastERA  as  a  flight  test  parameter 
identification  tool  will  be  assessed  by  exploring  the  effects  of  window  length,  sample  rate,  and  noise 
level  on  the  algorithm.  Finally,  FastERA  is  used  to  identify  the  frequency  and  damping  modes  of  the 
accelerating  simulation  results  for  the  two  degree  of  freedom  model  of  the  F- 16  that  is  developed  in 
chapter  3.  The  results  are  compared  to  the  results  of  an  actual  F-16  providing  validation  for  both  the 
use  of  FastERA  for  accelerating  data,  and  the  accuracy  of  the  2  degree  of  freedom  model. 

The  results  that  were  obtained  from  the  comparison  of  frequency  and  damping  values  between 
flight  test  and  the  two  degree  of  freedom  model  of  the  F-16  convinced  the  Air  Force  Flight  Test  Center 
(AFFTC)  to  support  a  deeper  exploration  into  the  feasibility  of  the  accelerated  flutter  testing  approach. 
The  Air  Force  Test  Pilots  School,  located  at  the  AFFTC  flew  F-16  test  sorties  to  collect  both  stabilized 
and  accelerated  flight  test  data  for  a  configuration  consistent  with  3  dimensional  simulations  that  were 
being  prepared  at  the  University  of  Colorado.  The  details  of  the  flight  test  performed  by  the  Test  Pilots 
School  are  presented  in  chapter  5  and  the  3  dimensional  simulations  will  be  presented  in  chapter  6. 

The  third  issue,  namely  “the  need  to  anticipate  the  vibration  response  at  airspeeds  beyond  which 
the  aircraft  has  accelerated”,  will  be  addressed  in  Chapter  7.  An  analytical  study  of  a  one  degree  of 
freedom  system  and  a  two  degree  of  freedom  system  is  developed.  Estimators  are  developed  for  these 
systems,  which  are  shown  to  be  overly  conservative.  Finally,  a  symbolic  solver  is  used  to  get  an  exact 
solution. 

The  final  issue,  namely  “developing  formal  criteria  to  determine  when  the  pilot  needs  to  end  the 
acceleration  and  terminate  the  test  point”,  is  also  addressed  in  chapter  7.  Utilizing  the  solutions  for  the 
one  and  two  degree  of  freedom  systems  to  project  the  flutter  speed,  the  acceleration  can  be  decreased 
or  stopped  when  the  control  room  personnel  and  the  pilot  have  ample  time  to  react. 
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1.7  Contributions  of  this  Thesis 

The  contribution  of  this  thesis,  to  the  flight  test  community,  is  the  development  of  a  new  flight 
test  procedure  for  conducting  flutter  testing  for  advanced  military  aircraft.  This  new  flight  test 
procedure  will  expedite  flight  testing  by  providing  a  means  of  conducting  flutter  tests,  and  identifying 
frequency  and  damping  parameters  for  the  critical  modes  of  the  aircraft  during  an  accelerated  flight. 

The  proposed  flight  test  technique  raises  theoretical  and  algorithmic  issues  that  have  not  been 
addressed  previously.  Specifically,  the  effect  of  acceleration  on  the  flutter  parameters  of  the  aircraft, 
whether  parameters  can  be  identified  using  data  from  an  accelerating  aircraft,  and  the  real  time 
prediction  of  the  aircraft  damping  beyond  the  current  airspeed  needed  to  be  investigated. 

Therefore  the  contributions  of  this  thesis  to  the  academic  community  are: 

1)  The  quantification  of  the  effect  acceleration  has  on  an  aerodynamic  system.  It 
was  shown  that  this  effect,  although  real,  is  not  of  significant  magnitude  to  invalidate 
test  results. 

2)  The  parameter  identification  using  data  from  an  accelerating  aircraft  was 
accomplished  using  small  window  sizes  within  which  the  parameters  do  not  change 
significantly.  In  practice  the  acceleration  of  a  typical  aircraft  is  not  sufficient  to  alter 
the  results  of  the  parameter  estimation  within  the  small  window. 

3)  The  use  of  estimators  were  employed  in  order  to  predict,  or  extrapolate,  the 
damping  characteristics  at  higher  airspeeds,  this  thesis  develops  both  an  estimator 
and  an  exact  solution  for  a  simplified  aeroelastic  system  found  in  the  classical 
literature.  The  estimator  was  found  to  be  unacceptably  conservative,  and  the  exact 
solution  was  found  to  be  quite  cumbersome  even  for  the  simplified  system. 


However,  it  is  shown,  at  least  in  principle,  that  the  accelerated  flutter  technique  is 
viable. 

A  preliminary  design  of  the  criterion  to  be  used  to  assure  that  the  testing  is 
halted  before  the  aircraft  crosses  the  flutter  boundary  is  developed  to  support  the  new 
flight  test  technique. 

Finally,  validation  is  performed,  primarily  to  validate  the  flight  test  technique. 
However,  this  validation  will,  as  a  consequence,  validate  the  simulation  technology. 
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2.0  Requirements  for  the  New  Flight  Test  Technique 


2.1  Safety  of  Test 

It  is  important  that  any  new  flight  test  methodologies  maintain,  or  increase  the  level  of  safety 
provided  by  the  methodologies  which  they  replace.  The  proposed  method  of  accomplishing  parameter 
identification  and  clearance  to  proceed  as  the  aircraft  accelerates  across  an  altitude  band  might  seem 
cavalier  at  first.  However,  there  are  safety  advantages  to  using  this  methodology.  It  is  a  truism  that 
most  flight  test  accidents  occur  when  the  aircraft  is  not  performing  a  test  point.  It  is  during  this  period 
when  the  pilot  and  the  control  room  are  more  relaxed  that  mistakes  can  more  easily  be  made. 

Therefore  there  is  a  safety  gain  because  several  stabilized  test  points  are  replaced  by  one  accelerating 
test  point. 

While  conducting  flutter  testing  the  data  sometimes  are  difficult  to  interpret  because  the 
resolution  of  the  stabilized  test  points.  For  example  if  a  downward  trend  in  damping  is  expected  at  a 
given  Mach  number,  say  0.84  Mach,  but  in  reality  the  trend  starts  0.83  Mach.  If  you  have  a  data  point 
at  0.80  Mach  and  0.85  Mach  it  will  not  be  clear  whether  the  break  occurred  at  0.83  Mach  or  at  0.84 
Mach  and  was  much  more  severe  than  anticipated.  This  confusion  is  compounded  by  noise  in  the 
damping  measurement.  Instead  of  points  being  exactly  identified  there  is  some  confidence  region 
around  the  points.  In  practice  when  these  kinds  of  uncertainties  suggest  a  potential  hazard  more  test 
points  will  be  added  to  better  characterize  the  region.  With  the  proposed  methodology  the  parameter 
identification  would  be  continuous  so  there  would  be  another  safety  gain  by  using  this  methodology. 
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The  simulation  being  used  in  this  research  is  higher  fidelity  than  has  previously  been  used  to 
support  flutter  flight  tests.  These  high  fidelity  simulations  alone  result  in  higher  confidence  in  the 
predictions  which  subsequently  results  in  a  higher  level  of  safety  during  test. 


2.2  Noise  Tolerance 

There  are  multiple  sources  of  noise  that  affect  flight  test  data.  First,  there  is  atmospheric 
turbulence  which  is  not  accounted  for  in  most  data  reduction  algorithms.  Second,  between  the  engines 
and  the  airflow  an  airplane  represents  a  fairly  harsh  vibration  environment.  Third,  the  input  forces 
cannot  be  measured  directly.  Instead  the  deflection  of  the  control  surface  is  measured.  Fourth,  all  of 
the  usual  measurement  errors  apply  to  aircraft  instrumentation  as  well. 

Because  flight  test  data  are  inherently  noisy,  it  is  critical  that  any  data  analysis  method  that  is  used 
is  tolerant  of  noise.  In  other  word  the  data  analysis  method  must  be  capable  of  identifying  frequencies 
and  damping  values  even  when  there  is  a  high  percentage  of  noise  in  the  data. 


2.3  Sample  Rate  Effects 

One  method  of  making  algorithms  more  tolerant  of  noise  is  to  collect  more  data.  A  straight 
forward  way  to  collect  more  data  is  to  increase  the  rate  at  which  the  data  are  sampled.  When  an 
algorithm  is  used  that  depends  on  Fourier  transforms  the  increased  sample  rate  results  in  decreased 
frequency  resolution.  This  may  in  turn  adversely  affect  the  algorithms  ability  to  correctly  identify  the 
frequency  and  damping  values. 

There  is  also  a  practical  limit  to  the  rate  at  which  the  data  can  be  sampled  and  telemetered  to  the 
ground.  The  government  continues  to  sell  radio  frequencies  to  private  industry.  With  the  proliferation 
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of  wireless  technology,  bandwidth  that  has  traditionally  been  used  to  support  flight  test  telemetry  is 
lost.  For  this  reason,  sample  rates  for  flutter  testing  will  likely  remain  at  the  current  100  -  800  samples 
per  second. 

Therefore  it  is  a  requirement  that  the  parameter  identification  algorithm  used  can  identify 
frequency  and  damping  values  from  data  at  these  sample  rates. 


2.4  Real  Time  Requirement 

It  is  a  requirement  that  frequencies  and  damping  values  be  calculated  real  time  as  well  as 
predictions  and  knock-it-off  calculations  be  performed  real  time.  Otherwise,  accelerating  the  aircraft 
as  the  flutter  analysis  is  performed  would  pose  too  great  of  a  safety  risk. 
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3.0  Effects  of  Acceleration  on  Structural  Parameters 

It  is  important  when  considering  the  use  of  flight  test  data  from  an  accelerating  aircraft  to  verify 
flutter  parameters  remain  essentially  the  same,  as  they  would  be  for  an  aircraft  in  stabilized  flight  at  the 
same  conditions.  For  this  reason  this  chapter  examines  the  effects  of  the  acceleration  on  the  flutter 
parameters.  This  issue  will  be  investigated  by  an  analytical  study.  The  analytical  model  depends  on 
the  equations  of  motion  that  will  be  presented  for  a  two  degree  of  freedom  “typical  wing  section”. 

Next  the  equations  of  motion  will  be  modified  to  account  for  the  effects  of  accelerating  the  frame  of 
reference,  which  is  attached  to  the  fuselage  of  the  aircraft.  Finally,  these  equations  will  be  applied  to  a 
two  dimensional  representation  of  an  F-16  wing,  and  it  will  be  shown  that  for  an  F-16  the  effect  of 
acceleration  will  be  negligible.  The  primary  purposes  of  this  chapter  are,  to  quantify  the  effects  of 
acceleration  on  the  structural  parameters,  and  to  demonstrate  that  those  effects  are  not  significant  for  a 
realistic  aircraft  conditions. 

3.1  Discrete  Equations  of  Motion 

The  equations  of  motion  for  a  typical  wing  section  [2]  (see  Fig.3.1)  in  the  absence  of  accelerations 
can  be  written  as  follows. 


m  h+  S0  0+  khh  =  Qh 


(3.1) 
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S0  h+ 10  6+  Ke6  -  Qe  (3.2) 


Where: 

m  total  mass  per  unit  span  of  the  wing 
h  bending  degree  of  freedom 

6  torsional  degree  of  freedom 

Se  m(xQ  +  )  static  moment  of  inertia 

IB  polar  moment  of  inertia 

Kh  bending  stiffness  coefficient 

Ke  torsional  stiffness  coefficient 

Qh  resultant  aerodynamic  force 

Qe  resultant  aerodynamic  moment 

G  center  of  gravity 

C  center  of  elasticity 


Figure  3.1:  A  typical  wing  section: 
elastic  center  (C),  center  of  gravity 
(G),  fuselage  (F),  angle  of  attack  ( 6 ). 
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Because  the  reference  frame  is  attached  to  the  aircraft  for  equations  (3.1)  and  (3.2)  inertial  terms 
must  be  added  whenever  the  aircraft  undergoes  acceleration  at  the  center  of  gravity.  When  the  aircraft 
is  accelerating,  but  6  remains  small  then  the  dynamic  equations  of  motion  become: 


mh+Se0+khh  =  Ql  (3.3) 

Seh+Ied+K*eO  =  Q*e  (3.4) 

Where  the  inertial  terms  that  result  from  the  acceleration  of  the  aircraft  are  given  by  the  equations 
below. 


K*e=Ke-  my Fx  (4  -  4 )  -  my  Fy  (y°c-y°c)  (3.5) 

Ql  =  Qh~mrFy  (3-6) 

Qe = Qe + mY Fx(y  g  ~  y c)  ~  my ^0% — xc)  (3.7) 


Where 

F  Denotes  a  fixed  point  on  the  fuselage 

y F  Acceleration  in  the  x  (horizontal)  direction  (note:  positive  denotes  deceleration  of  the  wing 
section) 

yF  Acceleration  in  the  y  (vertical)  direction 

ry 
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The  superscript  0  denotes  the  value  at  time  t  -  0 . 

As  can  be  seen  from  Equations  (3.3)  and  (3.4)  above,  accelerating  the  wing  section  changes  both 
of  the  aerodynamic  forces  and  the  torsional  stiffness.  These  changes  are  simply  the  result  of  inertial 
forces  due  to  the  fact  that  the  equations  represent  an  accelerating  system  and  a  frame  of  reference  that 
is  fixed  to  the  aircraft. 

The  torsional  force  changes  because  the  center  of  elasticity  and  the  center  of  gravity  do  not 
coincide.  As  a  result  the  inertial  force  acts  on  this  moment  arm  to  produce  a  torsional  force. 

If  the  airfoil  rotates,  as  shown  in  figure  3.2,  from  its  original  position  the  moment  arm  changes  in 
length  (for  small  angles  of  rotation  the  moment  arm  changes  by  (x°  -x°c)a)  resulting  in  a  torsional  force 

that  is  proportional  to  the  rotational  displacement.  It  is  this  term,  myF  -x£)  >  that  is  added  to  the 
torsional  stiffness  term,  Ko . 


,Wfx 

1}1Yfx 


-4 - 

Figure  3.2:  Torsional  stiffness  is  affected  by 
acceleration  because  the  inertial  force  my F  acts 
on  the  moment  arm  (*°  -x°c)a  (for  small  angles 
a). 
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The  result,  for  an  aircraft  that  has  a  forward  acceleration,  the  stiffness  increases  if  the  center  of 
gravity  is  aft  of  the  center  of  elasticity  due  to  the  term  myF  (x°  -x° )  .  Similarly,  if  the  aircraft  has  an 
upward  acceleration  component  the  rotational  stiffness  will  be  increased  by  center  of  gravity  is  below 
the  center  of  elasticity  due  to  the  term  myF  (y^G  —  y°c  ) .  However,  for  level  accelerated  flight  the 

acceleration  rate  yF  is  small  (perhaps  0.01  Mach/sec),  and  the  term  (x°  -  x°c  )  will  also  be  small 
because  x°G  and  x°  are  normally  located  within  one  quarter  chord  of  each  other.  Therefore,  for  most 
if  not  all  aircraft  the  quantity  my F  (x°  -x°  )  is  negligible.  Also  for  accelerated  level  flight  yF  =  0 . 

Therefore  it  can  be  reasonably  argued  that  K*0  =  Ke .  All  other  coefficients  of  the  left  hand-side  of 

Equations  (3.6)  and  (3.7)  are  unaffected  by  the  acceleration.  Therefore  it  can  be  concluded  the 
aeroelastic  parameters  (i.e.  frequency  and  damping)  for  level  accelerated  flight  are  the  same  as  for 
accelerated  flight  conditions. 

So  to  conclude  there  are  some  minor  effects  that  acceleration  has  on  the  underlying  parameters  of 
the  system.  However,  in  practice  these  effects  are  so  small  they  can  be  ignored  without  consequence. 
Section  3.2  will  use  an  airfoil  model  that  represents  an  F-16  to  demonstrate  the  acceleration  effects  can 
be  ignored. 


3.2  A  Typical  (2-D)  Wing  Section  Model  for  the  F-16 

This  model  was  developed  to  match  the  dynamic  characteristics  of  the  F-16.  This  model  will  be 
used  to  demonstrate  the  effect  of  acceleration  on  the  flutter  parameters.  In  chapter  4  simulations  will 
be  performed  using  this  model  in  order  to  validate  the  use  of  ERA  and  the  Accelerated  Flutter  Test 
technique. 
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Starting  from  a  detailed  finite  element  structural  model  of  a  “clean”  right  wing  of  the  F-16  Block 
40  equipped  with  a  missile  launching  system  at  its  tip  (Fig.  3.3).  The  objective  is  to  construct  a  two- 
degree  of  freedom  wing  section  model  that  is  equivalent  to  the  three-dimensional  wing  in  the  following 
sense 

1 .  It  reproduces  the  first  bending  and  torsion  frequencies,  which  are  predicted  by  the  three- 
dimensional  finite  element  model  of  the  wing  to  be  4.76  Hz  and  7.43  Hz,  respectively. 

2.  It  reproduces  the  same  bending  and  torsion  modal  masses. 

3.  It  reproduces  the  same  vertical  displacement  at  the  leading  edge  of  the  cross  section  located  at 
68%  of  the  distance  from  the  root  to  the  tip  of  the  three-dimensional  wing,  for  both  bending 
and  torsion  mode  shapes,  when  these  are  normalized  to  a  unit  rotation  of  the  cross  section. 

When  the  wing  has  a  large  aspect  ratio  and  a  small  sweep  angle,  it  is  commonly  suggested  to 
choose  the  mechanical  properties  of  the  typical  wing  section  as  to  match  those  of  the  cross  section 
located  at  70%  to  75%  of  the  distance  from  the  root  to  the  tip  of  the  three-dimensional  wing  (for 
example,  see  [2]  and  references  6-1  and  6-2  therein).  Nevertheless,  the  cross-section  located  at  68%  of 
the  distance  from  the  root  to  the  tip  (see  Fig.  3.3)  was  chosen,  as  stated  in  criterion  3  above,  because 
the  F-16  wing  is  strongly  tapered  and  is  rather  soft  towards  its  tip. 
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1 :  cross-section  chosen  for  the  structural  properties 
2:  cross-section  chosen  for  the  aerodynamic  properties 

Figure  3.3:  Three-dimensional  detailed  finite  element  model  of  an  F-16  wing 


Besides  the  shape  of  the  airfoil,  six  parameters  define  the  sought-after  typical  wing  section,  These 
are  denoted  collectively  by,  Pr  s  where  the  subscript  tws  stands  for  typical  wing  section,  and  listed 
below 


P,,s  =  {mJo>x  o-TV; -xc),Kh,K0} 


(3.8) 


The  three  criteria  stated  above  for  establishing  the  equivalence  between  the  typical  wing  section 
and  the  three-dimensional  wing  can  be  formulated  as  follows 


2ttx  A.lGrad  /  5 

^,,2(Tv) 

2k  x  lA3rad  /  5 

1.375xl06^g-m2 

Mnvs,2(Pnvs) 

2.523 Kg -in¬ 

is. 6\m 

JtAP,s)_ 

1.017/;? 

(3.9) 


36 


Where  D.  denotes  a  frequency,  the  subscripts  1  and  2  refer  to  the  bending  and  torsion  modes, 
respectively,  //  denotes  a  modal  mass,  the  superscript  le  designates  the  leading  edge, 

h'e  =  h  +  ( x*  -  xc  ){Q  =  1)  is  the  vertical  displacement  at  the  leading  edge  of  the  typical  wing  section  and 
all  components  of  the  right  hand-side  of  Equation  (3.9)  are  obtained  from  the  detailed  finite  element 
model  of  the  three-dimensional  wing.  Note  that  the  fact  that  hl^sX  <£C  h^vs  2  suggests  that  despite  the 

bending/torsion  coupling,  the  first  mode  of  the  wing  is  dominated  by  bending,  while  the  second  one  is 
dominated  by  torsion. 

The  constraints  (3.9)  lead  to  a  nonlinear  system  of  equations  with  six  unknowns,  which  is  solved 
by  the  Nelder-Mead  simplex  method  (function  fminf  in  Matlab)  to  obtain 

m  =  2.05xl03Xg 
I0=2.53x\O3Kg-m2 

XG  =1.1 26m  (upstream  of  section  at  68%) 

(xG  -  Xc  )  =  0 . 0642m  (G  behind  C) 

Kh  =  2.046  xl06rV7m 
Kd  =  5.468  xl06iV7m 

Using  the  above  numbers,  the  modified  torsional  stiffness  (see  Equation  (3.6))  for  an  acceleration 
as  high  as  0.05  Mach/second  is  Ke  =  5.471  x  106  Nm.  Hence,  as  expected,  KQ  =  K& ,  which  supports 

the  conclusion  made  above  regarding  the  negligible  effect  of  level  flight  acceleration  on  the  aeroelastic 
parameters  of  a  system. 
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4.0  System  Identification  using  ERA 

Depending  on  several  factors,  among  which  is  the  flow  regime,  an  aeroelastic  system  can  behave 
linearly  or  non-linearly.  This  raises  a  first  albeit  minor  concern  as  to  the  applicability  of  conventional 
signal  processing  techniques  to  the  identification  of  the  aeroelastic  parameters  of  an  aircraft,  or  a 
typical  wing  section  particularly  in  the  transonic  regime.  Furthermore,  an  accelerated  aeroelastic 
system  is  also  a  time-varying  system.  This  is  essentially  because  the  mass,  damping  and  stiffness 

properties  of  the  “wet”  structure  are  influenced  by  the  free  stream  velocity  Vm  of  this  structure,  and 

varies  with  time  during  an  accelerated  flight.  This  raises  a  second  concern  as  to  the  continuous 
parametric  identification  of  an  accelerating  aircraft  using  existing  signal  processing  analysis  which,  in 
theory,  are  limited  to  time  independent  systems.  Both  concerns  can  be  addressed  by  the  windowing 
approach  described  below,  which  is  in  principle  applicable  to  several  identification  methods. 

First  note  that: 

•  some  methods  are  capable  of  identifying  the  frequency  and  damping  coefficient  of  the  lowest 
mode  of  a  structure  using  as  few  as  2  cycles  of  the  response  of  this  structure  to  an  input  perturbation 

•  for  a  mode  at  10  Hz,  2  cycles  corresponds  to  a  time  interval  of  0.2  second. 

•  at  an  altitude  of  10,000  ft  a  level  flight  acceleration  of  0.05  Mach/second  corresponds  to  a 


horizontal  acceleration  of  1 .6  g’s.  Such  acceleration  is  beyond  the  reach  of  most  if  not  all  aircraft. 
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•  during  a  time-interval  of  0.2  second,  the  speed  of  an  aircraft  accelerating  at  0.05  Mach/second 
varies  by  0.01  Mach. 

•  a  1%  variation  of  a  Mach  number  has  a  negligible  effect  on  the  frequency  and  damping 
coefficient  of  a  wet  structure  cruising  at  . 

From  the  above  observations  it  can  be  concluded  that  within  a  time-window  of  the  order  of  0.2 
second,  an  aeroelastic  system  with  a  first  wet  mode  at  1 0  Hz  can  be  considered,  for  all  practical 
purposes,  as  a  time-invariant  system.  Furthermore,  within  a  window  of  that  size  the  aeroelastic 
response  of  such  a  system  can  be  assumed  to  be  linear,  as  long  as  the  structure  itself  behaves  linearly, 
which  is  usually  the  case  for  an  aircraft  excited  by  an  input  to  a  control  surface.  Hence  it  can  also  be 
concluded,  at  least  in  principle,  that  it  should  be  possible  to  expand  the  flutter  envelope  of  an  aircraft 
using  its  continuous  vibration  response  to  input  perturbations  during  an  accelerated  flight  by  applying 
the  following  simple  procedure 

1)  locate  the  time  instances  at  which  the  Mach  numbers  of  interest  are  reached  by  the  aircraft. 

2)  for  each  Mach  number  of  interest,  define  a  time-window  of  width  equal  approximately  to  0.2 
second  or  two  cycles  of  the  expected  lowest  frequency. 

3)  within  each  time-window,  apply  a  suitable  signal  processing  method  to  identify  the  frequency 
and  damping  coefficient  of  the  wet  structure  and  associate  these  with  the  Mach  number  for  which  this 
window  is  defined. 

The  window-based  identification  approach  summarized  above  can  be  performed  either  real-time 
or  off-line  as  a  post-processing  signal  analysis  tool.  In  this  work,  the  Eigenvalue  Realization 
Algorithm  (ERA)  [3]  is  chosen  to  perform  the  parameter  identification. 
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The  ERA  is  an  identification  method  for  linear  time-invariant  systems.  It  implicitly  assumes  that 
the  dynamic  response  of  the  given  system  is  sampled  at  a  constant  rate.  It  can  handle  multi-degree  of 
freedom  systems,  and  is  less  sensitive  to  noise  than  the  logarithmic  decay  method.  In  order  to  keep 
this  chapter  self-contained,  an  overview  of  ERA  will  be  presented  along  with  its  fast  implementation 
FastERA  [4,5]  which  is  suitable  for  real-time  processing. 


4.1  Overview  of  ERA 

In  a  linear  context,  damped  structural  vibrations  are  governed  by  the  following  equations  of 
dynamic  equilibrium. 


••  • 

Mq+Dq+Kq  =  f(t )  (4.1) 

where  M ,  D ,  K  are  respectively  the  mass,  damping,  and  stiffness  matrices,  and  where  / ( t )  denotes 
the  vector  of  external  forces.  The  above  equations  can  be  re-written  in  state-space  form  as  follows. 


x  =  kx  +  Bf(t ) 


(4.2) 


q  =  Cx  (4.3) 


where 


(4.4) 


K  = 


0  / 
-M~XK  -M~XD 


B 


0 

M~x 


C  =  [I  0] 
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Assuming  that  initially  the  structural  system  is  at  rest,  and  that  it  is  excited  at  t  =  0  by  an  impulse 
load  f(t)  =  u0S(t)  where  S(t)  denotes  the  Dirac  delta  function,  the  solution  of  Equations  (4.2)  and 
(4.3)  is: 

q(t )  =  CeKlBu0  (4.5) 

Note  that  Equation  (4.5)  is  also  the  solution  of  Equations  (4.2)  and  (4.3)  for  f(t)  =  0  and  an 
initial  velocity  equal  to  M~lu0 . 

For  an  arbitrary  excitation,  / (/)  can  be  represented  by  a  series  of  impulses  at  discrete  time- 

oc 

stations  t{-  that  is  /(/)  =  ^  u.S(t  - 1.)  ,  in  which  case  the  solution  of  Equations  (4.2)  and  (4.3)  at  time 

»=i 

tk  is  given  by 

qk  =  ?(/,)  =  £  Cev>  =  £  Ce(k~')£JK  Bu,  (4,6) 

1=0  /= 0 

where  a  constant  sampling  rate  V  At  is  assumed  so  that  t k  - 1.  =(k-  i)At . 

Let 

A  =  eAtK  (4.7) 

Using  the  above  definition  of  the  matrix  A  ,  Equation  (4.6)  can  be  re-written  as  follows. 
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(4.8) 


qk  =  ^  CAk~‘Bui  =  YjCAk~i-l(AB)ul+CBuk  =  Yj{CAk~i-xBui)  +  Guk 

i= 0  i=0  /= 0 


where 


G  =  CB  and  B  =  ^5  (4.9) 

Hence  qk  is  the  classical  solution  of  the  state-space  realization  problem  for  the  following  linear 
sampled  data  system. 


xk+\  =  Axk  (4.10) 

<Hk=CXk+Gllk k  (4.11) 

If  ndoj-  denotes  the  total  number  of  degrees  of  freedom  of  the  computational  model  described  by 

equation  (4.1),  but  only  n^f  <  ndof  degrees  of  freedom  are  measured,  then  the  matrix  C  can  be  written 
as: 


C  =  [L  0]  (4.12) 

where  L  is  an  n™  x  ndof  Boolean  matrix.  Hence,  in  general,  the  number  of  state-space  variables  Xk 
that  are  considered  is  nx  =  2 ndof  ,  and  the  matrix  A  defined  in  Equation  (4.7)  is  an  nx  x  nx  matrix. 


In  general,  the  discrete  convolution  sum  (4.8)  is  expressed  as: 
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(4.13) 


where  the  matrices 


Mk  {cAm~]  B 

are  known  as  the  Markov  parameters.  Hence,  qk  can  be  directly  related  to  the  load  components  u.  via 
the  Markov  parameters  A/ . 

Let  zs  and  —  as  ±  icos  denote  the  complex  eigenmodes  and  eigenfrequencies  of  the  structure 
defined  by 


9m  =  0 
9m  >  0 


(4.14) 


X*Mzs  +  AsDzs  +  Kzs  —  0  (4.15) 

From  Equation  (4.4),  it  follows  that  the  A  are  also  the  eigenvalues  of  K ,  and  the  eigenvectors  of 
this  matrix  are 


^5  = 


(4.16) 


The  ERA  exploits  the  results  summarized  above  as  follows.  First,  it  constructs  the  Markov 
parameters  Mk  using  the  input  and  output  data,  namely,  the  impulse  loads  and  measured 


displacements.  Then,  it  extracts  the  matrix  A  from  the  Markov  parameters  M k  using,  for  example, 
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the  FastERA  algorithm  described  in  Section  4.2.  Finally,  the  ERA  computes  the  eigenvalues  <r.  and 
eigenvectors  as  of  A  which  satisfy 


Aas=crsas  (4.17) 


From  Equation  (4.7),  it  follows  that  the  sought-after  complex  modes  ,  damping  ratios  £ ,  and 
frequencies  to  of  the  given  structural  system  are  given  by 


zs  =  Cas 


4ARe{lnM 


0)c 


=^7Im{ln(°v)} 


(4.18) 


4.2  The  FastERA  Algorithm 

This  work  uses  the  FastERA  implementation  or  the  ERA  to  perform  all  identifications.  For  this 
reason,  FastERA  is  reviewed  in  this  section  and  some  key  variables  are  introduced  that  appear  in  the 
remainder  of  this  work. 


The  FastERA  method  is  based  on  the  analysis  of  the  following  Hankel  matrix  that  is  defined  for  a 
data  set  (qk ,  uk )  sampled  at  N  points  in  time 


Mx  M2  •••  Md 

m2  m3  ... 


Mq  Mq+X 


M 


q+d+ 1 


(4.19) 
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and  where  q  and  d  are  such  that  N  =  q  +  d  .  From  the  definition  (4. 14)  of  the  Markov  parameters, 
it  follows  that  Hqd  can  be  decomposed  as  follows 

Hqd=VqWd  (4.20) 


where 


c 


CAq~l 


Wd=[B  AB  •••  Ad~lB 


(4.21) 


The  qn™f  x  qndoj-  matrix  Vq  determines  the  observability  of  the  system  whereas  Wd  is  related 
to  its  controllability. 


The  FastERA  method  starts  with  the  construction  of  the  qi ld^  x  qndoj-  “square  data  matrix” 


Jq=HqdHTqd=VqWdWjv;  (4.22) 

For  a  given  set  of  inputs  uk  and  outputs  qk  ,  the  choice  of  state-space  variables  X  is  not  unique. 

Indeed,  the  input/output  relation  is  not  affected  by  the  change  of  variable  X  =  T  X ,  which  transforms 
the  realization  system,  (4.10)  and  (4.1 1)  into 

xk+\  =T~'ATxk  +  TlBuk 


qk  -CTxk  +  Guk 


(4.23) 
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The  controllability  factors  associated  with  x  are 

_  T 

wd  Wd  =  T-xwdwJrT  (4.24) 

T 

In  particular,  choosing  T  =  (WdWj)V2  leads  to  WdWd  -  I .  Hence,  without  any  loss  of 

generality,  the  FastERA  algorithm  assumes  that  the  state-space  variables  are  chosen  so  that  Equation 
(4.22)  simplifies  to 

(4-25) 

From  Equation  (4.21),  it  follows  that  the  first  (q  —  rows  of  Vq  are  given  by 


(4.26) 


and  the  last  ( q  —  1  rows  of  that  matrix  are  given  by 

'  CA  " 

vm=  CA’  =vmA  ,4.27, 

CAqA 
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In  general,  the  number  of  degrees  of  freedom,  ndoj  associated  with  a  system  to  be  identified  is 
unknown  a  priori.  Indeed,  tldoj  is  infinite  for  any  continuous  system.  In  practice  the  target  values  for 

ndof  and  Ae  corresponding  number  of  state-space  variables  nx  =  2 ndoj-  are  dictated  by  the 
complexity  of  the  model  to  be  realized.  Once  the  user  sets  these  values,  the  FastERA  algorithm 
computes  the  nx  largest  eigenvalues  K.  of  J q  and  their  corresponding  eigenvectors  P-  in  order  to 
build  the  square  root  factorization 

Jq  =  Pdiag(Kj)PT  =  {Pdiag(K)2){diag{KX!2)PT)  =  Vq  Vq  (4.28) 


where  P  =  [/>,  •  •  •  p.  •  •  •]  and  Vq  =  Pdiag(ic J'2 ) .  From  Eqs.  (4.27)  applied  to  V  ,  it  follows  that 


A  = 


(  ~(1)V  ~(2) 
Vq  Vq 
v  J 


(4.29) 


(~0)) 

♦ 

where 

Vq 

is  the  pseudo-inverse  of 

V, 

K  J 

Hence,  the  extraction  of  the  nx  largest  eigenpairs  of  the  x  qn ^  data  matrix  J q  is  the 

most  computationally  significant  step  of  the  FastERA  algorithm.  Keeping  q  relatively  small,  say 
q~d!  5  allows  FastERA  to  operate  in  real  time.  For  further  details  on  this  identification  method, 


refer  to  [4], 


Frequency  Relative  Error  {%) 
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Next,  the  window-by-window  application  of  the  ERA  is  validated  for  the  identification  of  the 
parameters  of  an  accelerating  aircraft  by  a  series  of  numerical  simulations  designed  for  an  F-16  fighter 
configuration  for  which  flight  test  data  are  available. 


4.3  ERA  Sample  Rate  Effects 

It  is  not  only  expensive  to  telemeter  large  amounts  of  data  to  the  control  room,  but  with  multiple 
programs  competing  for  available  bandwidth,  and  the  decrease  in  available  band  width  due  to 
frequencies  being  sold  for  commercial  applications,  sometimes  it  is  not  even  possible.  Because  of  the 
need  to  minimize  the  amount  of  data  that  must  be  telemetered  to  the  control  room  a  study  was 
conducted  to  determine  the  effects  and  limitations  of  reducing  the  sample  rate.  Fig.  4.1  shows  the 
effect  that  various  sample  rates  have  on  ERA’s  ability  to  accurately  identify  the  frequency  of  the 
torsional  bending  modes.  This  figure  shows  that  for  sample  rates  below  about  20  samples  per  cycle 
the  algorithm  cannot  accurately  identify  the  damping  ratio.  For  the  frequency  (see  Fig.  (4.1))  the 
algorithm  cannot  accurately  identify  below  about  8  samples  per  cycle. 


Effect  of  Sampling  Rate  on  Identified  Frequency  (Modes  1  and  2) 


10  10  10 
Number  of  Sampling  Points  per  Cycle 


Effect  of  Sampling  Rate  on  Identified  Damping  Ratio  (Modes  1  and  2) 


Figure  4.1 :  Effects  of  sampling  rate  on  parameters  from  noise  free  data 


48 


This  implies  that  for  clean  data  (data  with  no  noise)  the  minimum  acceptable  sample  rate  for  a  ten 
hertz  mode  would  be  about  200  cycles  per  second.  This  is  marginally  within  the  range  of  sample  rates 
normally  used  structural  envelope  expansion  flight  testing.  However,  this  is  for  clean  data,  and  may 
not  work  for  data  with  noise. 


4.4  ERA  Noise  Tolerance  Study 

Realizing  that  the  ERA  algorithm  is  based  on  Fast  Fourier  Transforms  and  all  the  problems  with 
using  small  block  sizes  apply,  namely  frequency  resolution  and  noise  tolerance,  a  noise  tolerance  study 
was  conducted  in  order  determine  the  usefulness  of  ERA  for  determining  frequency  and  damping 
characteristics  of  noise  ridden  data. 

Noise  was  generated  for  this  study  using  the  random  function  in  Matlab.  The  random  function 
produces  data  that  have  a  level,  but  not  constant  frequency  content.  Noise  levels  of  1%  and  5%  were 
used.  The  1%  noise  level  means  the  highest  possible  value  of  the  random  function  is  set  to  1%  of  the 
highest  excitation  level  in  the  signal.  5%  noise  is  defined  in  same  fashion. 

The  diagram  in  Fig.  (4.2)  shows  the  scatter  of  frequency  and  damping  verses  the  number  of 
samples  per  cycle.  For  the  damping  ratio,  the  scatter  stabilizes  at  about  !4%  error  at  about  1,000 
samples  per  cycle.  For  the  frequency,  the  scatter  stabilizes  at  about  %%  scatter,  and  that  again  is  at 
about  the  1,000  sample  per  cycle.  This  level  of  uncertainty  might  be  acceptable,  however  the  sample 
rate  would  have  to  be  5000  samples  per  second.  It  is  possible,  but  not  practical,  to  sample  and 
telemeter  flight  test  data  real  time  at  that  sample  rate. 
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Number  of  Sampling  Points  per  Cycle  Number  of  Sampling  Points  per  Cycle 

Figure  4.2:  Effects  of  sample  rate  on  parameters  from  data  with  5%  noise 


4.5  ERA  Window  Size/Acceleration  Rate  Study 


As  a  result  of  this  inability  to  satisfactorily  identify  modes  with  a  window  of  0.2  seconds  and 
noisy  data,  it  was  determined  that  the  window  size  should  be  increased.  But  the  block  size  is  limited 
by  the  acceleration  rate.  If  the  Mach  number  changes  too  much  in  the  window  then  the  algorithm  will 
fail.  So  the  above  noise  analysis  was  conducted  with  a  block  size  of  1 .0  second  and  the  acceleration 
rate  was  reduced  to  0.01  Mach  per  second. 


Mode  2:  Effect  of  Window  Length  on  Identified  Frequency  Mode  2:  Effect  of  Window  Length  on  Identified  Damping  Ratio 


Figure  4.3:  Effects  of  acceleration  rate  and  window  size  on  parameters  from  data  with  5%  noise 
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The  results  shown  in  Fig.  4.3  indicate  the  one-second  block  is  sufficient  to  identify  the  parameters 
even  in  the  presence  of  5%  noise. 


4.6  Validation 

There  are  two  things  that  need  to  be  validated  here.  They  are  the  two  issues  that  have  been 
examined  thus  far  in  this  thesis.  First,  that  acceleration  does  not  significantly  affect  the  aeroelastic 
parameters.  Second  that  the  windowing  is  an  effective  way  to  identify  the  aeroelastic  parameters. 

Both  of  these  validations  can  be  accomplished  simultaneously  by  comparing  simulation  data  from 
stabilized  conditions  to  simulation  data  from  accelerated  conditions.  Some  critics  might  not  be 
satisfied  with  this  approach  because  of  the  potential  that  both  simulations  are  flawed.  In  order  to 
address  this  concern  the  data  from  the  simulations  will  be  compared  to  flight  test  data  provided  by  the 
Air  Force  Flight  Test  Center.  Section  4.6.1  will  outline  the  details  of  the  modeling  process  for  the 
simulations  presented  in  this  chapter,  and  the  simulation  technology  that  is  used  throughout  this 
document.  Details  pertaining  to  how  the  simulated  structure  was  excited,  and  how  exploiting  the 
moving  mesh  algorithm  enforced  the  acceleration  of  the  model  will  also  be  discussed. 

There  are  minor  differences  in  configuration  between  the  simulations,  which  represent  a  clean 
wing  (i.e.  no  external  stores)  configuration,  and  the  flight  test  data  which  are  for  a  configuration  with 
minimal  stores.  Also  the  simulation  is  a  two-dimensional  representation  of  the  full  three-dimensional 
aircraft.  Even  so,  the  correlation  between  them  is  remarkable.  Figure  4.6  shows  that  predicted  wet 
frequencies  are  in  good  agreement  with  those  measured  in  flight  test.  The  predicted  damping  ratios  do 
not  agree  as  well  with  the  experimental  data,  even  though  the  trend  of  their  variation  with  the  Mach 
number  is  similar  to  that  of  the  flight  test  data.  As  a  flutter  flight  test  engineer  these  results  are  quite 
exciting  because  a  simple  model  is  providing  better  prediction  than  are  typically  available  to  support 
the  flight  test. 
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4.6. 1  Simulation  Technology 

Because  the  validation  is  dependent  on  simulation,  it  is  expedient  to  explain  the  technology  that 
was  used  to  perform  those  simulations.  Because  these  are  aeroelastic  simulations  it  is  implied  that 
there  is  some  representation  of  the  structure  and  some  representation  of  the  aerodynamics  (i.e.  the 
fluid).  And  the  technology  would  not  be  complete  without  discussing  the  three-field  approach  that  is 
used  to  manage  the  information  as  it  is  passed  between  the  structure  and  the  fluid  software  modules. 


4.6. 1 . 1  Structure  Model  Used  for  the  Simulations 

The  structure  that  was  used  in  these  simulations  is  the  typical  (2  degree  of  freedom)  wing  model 
representation  of  the  F-16  that  was  developed  in  chapter  3.  For  the  sake  of  brevity  this  model  will  not 
be  revisited  here. 


4.6.1. 2  Fluid  Model  Used  for  the  Simulations 

The  airfoil  of  the  typical  wing  section  was  chosen  as  that  of  the  cross  section  located  at  45%  of 
the  distance  between  the  root  and  tip  of  the  three-dimensional  wing  because: 

•  The  chord  of  the  airfoil  of  the  typical  wing  section  must  be  close  to  the  ratio  of  the  wetted  area 
and  the  wingspan. 

•  Because  the  wing  is  tapered,  most  of  the  lift  is  generated  by  the  section  of  the  wing  that  is  close 
to  the  root,  which  means  that  the  aerodynamic  center  of  the  wing  is  within  that  area. 
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The  flow  domain  was  discretized  by  18,000  vertices,  to  ensure  a  sufficient  resolution  for  shock 
capturing  in  the  region  close  to  the  sharp  leading  edge  (Fig.  4.4).  Because  the  purpose  of  the  typical 
section  is  to  represent  the  entire  wing,  each  aerodynamic  force  obtained  from  a  flow  computation  on 
this  two-dimensional  mesh  was  multiplied  by  the  span  of  the  wing. 
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4. 1.6.3  Discrete  Fluid  Equations  of  Motion 

To  adequately  describe  the  motion  of  the  typical  wing  section,  introduced  above,  requires 
evaluating  the  resulting  aerodynamic  force  and  moment  Qh  and  Qe  .  The  fluid  was  modeled  using  the 

Arbitrary  Lagrangian  Eulerian  (ALE)  form  of  the  Euler  equations  in  order  to  handle  the  vibrations  of 
the  airfoil  and  its  acceleration. 

Let  Q(/)  =  M2  be  the  flow  domain  surrounding  the  airfoil  of  the  typical  wing  section,  and  T(/) 
be  its  moving  boundary.  A  mapping  function  is  introduced  between  Q(/)  where  time  is  denoted  by  t 
and  a  grid  point’s  coordinates  by  X ,  and  a  reference  configuration  Q(7)  where  time  is  denoted  by  t 
and  the  grid  point’s  coordinates  by  £ ,  as  follows. 

X  =  x(^,z)\t-T  (4.30) 

The  ALE  conservative  form  of  the  Navier-Stokes  equations  describing  viscous  flows  on  dynamic 
meshes  can  be  written  as 


d(JW) 

— — -L  j7  (yy,x)  =  o  (4.3i) 

dt  h 


JTC  (W,  x)  -  f(W)  -xW  (4.32) 


•  Qx 

Where  a  dot  superscript  designates  a  time  derivative,  J  =  det( dx  /  d<%)  ,  X  = - 

dr 


is  the 


velocity  of  the  dynamic  fluid  mesh,  W  is  the  fluid  state  vector  (state  variables),  and  JFC  denotes  the 


ALE  convective  fluxes. 
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Semi-discretized  Equation  (4.31)  is  semi-discretized  on  a  triangulation  from  which  a  dual  mesh  is 
derived,  defined  by  time-dependent  control  volumes  or  cells  Ci(t )  .  Resolving  the  ALE  convective 
fluxes  by  a  suitable  Riemann  solver.  The  resulting  semi-discrete  equation  of  equilibrium  of  the  fluid  is 


^-(AlWl)  +  F,(W,X,X)  =  0  (4.33) 

at 


Where  Ai  =  dxQ ,  W{  denotes  the  average  value  of  W  over  the  cell  Ci(t) ,  F(  denotes 

CM 

the  semi-discrete  ALE  convective  flux,  W  is  the  vector  formed  by  the  collection  of  Wt ,  and  X  is  the 
vector  of  time-dependent  grid  point  positions.  Various  expressions  of  the  flux  approximation 

F.  ( W ,  X ,  X)  can  be  found  in  [6-9] 

Time  integrate  the  semi-discrete  equations  of  flow  motion  (4.33)  on  dynamic  meshes  using  the 
generalized  second-order  implicit  backward  difference  scheme  developed  in  [10,1 1],  This  scheme 
satisfies  the  second-order  discrete  geometric  conservation  law,  and  retains  second  order  time-accuracy 
on  moving  grids  [10,1 1].  It  can  be  written  in  compact  form  as  follows. 


(4.34) 


«„+i4 


n+1W”+1  +anA"Win+a„_lA1n-lWin-l+AtnFi(Wn+l,XnJel,X)  =  0 


.  n+1 

Where  —  A.(X‘ n)  >  an^  ,  and  Xave  are  the  following  averaged  position  and  velocity  of 


the  fluid  moving  mesh. 


yn+\  _  n 
ave  ot'«+ 1 


(xn+xn+'^ 


5m  (X^+XX 

c  V  2  7 


(4.35) 
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a  j »  CCn  ,  and  a  are  time  dependent  because  a  variable  time  step  (\t(yi)  is  employed  and  are 
given  by 


A  f 


A 


(4.36) 


«„-l  = 


_  (T2 


(4.37) 


a„=-l-C  (4.38) 


a 


n-\ 


=  l  +  2£ 

■  1+C 


(4.39) 


It  remains  to  specify  how  Xn ,  the  vector  position  of  the  fluid  grid  points,  is  updated  at  each  time- 
station  tn  . 

For  this  purpose,  first  note  that  for  a  typical  wing  section  problem,  the  simplest  strategy  for 
updating  the  position  of  the  fluid  mesh  is  to  move  it  rigidly  with  the  airfoil.  Then  note  that  the  motion 
of  the  airfoil  is  completely  determined  by  the  motion  of  the  fuselage  point  F  ,  and  the  vibrations  of  the 
bending  and  torsional  degrees  of  freedom,  h  and  6 .  It  follows  that  at  each  time-station  tn  ,  the 
position  of  any  fluid  grid  point  N  is  given  by: 


xnN=xnF+  ( x°N  - x°F)cos 6n  - (y°N  - y°F) sin 6n 


(4.40) 
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ynN=ynF+hn+  ( y°N  - yF ) cos 6n  - (x°N  - xF) sin 6n  (4.41) 

Where  hn  and  6n  are  determined  by  solving  the  coupled  fluid/structure  equations  of  motion 
(3.3),  (3.4),  and  (4.33),  and  the  instantaneous  position  ,  ynf  )  of  the  fuselage  point  F  is  deduced 
from  the  specified  acceleration  of  the  aircraft.  For  example,  for  a  constant  acceleration,  the 
instantaneous  position  of  the  point  F  is  given  by: 

XnF  =l^-tn2  +XF  (4.42) 

yt=^f+/F  («3) 

Hence,  the  acceleration  of  the  typical  wing  section  is  transmitted  to  the  fluid  by  Equations  (4.42), 

and  (4.43),  and  accounted  for  by  the  additional  convection  term  x  W  that  characterizes  the  ALE  form 
(4.31)  of  the  flow  equations. 

Remark.  At  the  first  glance,  the  reader  may  think  that  Equations  (4.42),  and  (4.43)  are  missing 
the  terms  Vm  tn  and  Vm  f  where  Vw  is  the  free-stream  velocity.  However,  these  terms  are  not 

missing.  They  are  automatically  taken  into  account  by  the  initial  conditions  of  the  CFD  simulation 
through  the  specified  free-stream  Mach  number  M ^  . 

The  coupled  fluid/structure  discrete  equations  of  motion  are  solved,  throughout  this  chapter,  by  the 
second-order  time-accurate  staggered  and  sub-iteration  free  algorithm  described  in  [12].  As  stated 
earlier,  this  staggered  algorithm  is  equipped  with  the  midpoint  rule  as  a  structural  time-integrator,  and 
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the  generalized  second-order  implicit  backward  difference  scheme  (4.34)  developed  in  [10,1 1]  as  a 
flow  time-integrator 


4.6. 1.4  The  Three  Field  Approach 

In  order  to  combine  the  computational  fluid  dynamics,  which  utilizes  a  coordinate  system,  fixed 
in  space  (Eulerian  coordinate  system),  and  the  finite  element  method,  which  utilizes  a  coordinate 
system,  fixed  to  the  material  (Lagrangian  coordinate  system),  it  is  critical  to  allow  the  computational 
fluid  dynamic  mesh  to  move  to  accommodate  the  motion  of  the  finite  element  mesh.  The  method  used 
to  allow  the  motion  of  the  fluid  mesh  is  to  formulate  the  mesh  as  a  psudo-structural  system.  The 
combination  of  the  computational  fluid  dynamics,  the  finite  elements,  and  the  psudo-structural  fluid 
mesh  system  is  known  as  the  three  field  formulation.  The  way  the  three  field  formulation  works  is 
through  message  passing  between  the  three  software  modules.  The  finite  element  module  sends 
displacement  information  to  the  mesh  motion  module  which  calculates  the  motion  of  each  grid  point  in 
the  fluid  mesh.  The  mesh  motion  is  then  sent  to  the  computational  fluid  dynamic  module.  The 
computational  fluid  dynamics  module  calculates  the  fluid  states  for  each  cell  of  the  fluid  grid.  The 
fluid  flux  related  to  the  motion  of  the  fluid  mesh  is  accounted  for  using  the  Arbitrary  Lagrangian 
Eulerian  (ALE)  formulation  described  in  reference  [13].  The  pressure  is  then  sent  to  the  finite  element 
module  and  the  cycle  is  complete.  The  timing  of  the  message  passing  is  important  to  maintain  the 
order  of  accuracy  of  the  simulation.  The  staggered  approach  described  in  reference  [14]  is  used  for 
this  reason.  The  equations  governing  each  algorithm  are  listed  below. 


d[jw ) 


dt 


+  JV« 


8x 


dt 


w 


=  0 


(4.44) 


Fluid  Equation 
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Structure  Equation 


d\ 
dt 2 


-div 


( 7 . 


)) 


(4.45) 


Mesh  Motion  Equation 


~  d2x 

p — t —  div 
H  dt2 


f 


w 


=  o 


j) 


(4.46) 


There  are  also  some  boundary  conditions  that  must  be  observed.  First  on  the  boundary  between 
the  fluid  and  the  structure  the  normal  pressure  of  the  fluid  has  to  equal  the  normal  stress  of  the 
structure  and  the  normal  velocity  of  the  fluid  has  to  equal  the  normal  velocity  of  the  structure.  Second 
the  normal  displacement  and  the  normal  velocity  of  the  fluid  mesh  and  structure  must  be  equal  along 
the  boundary  between  them. 


Fluid/Structure  Interface  Conditions 


—pn  =(7s-n  (4.47) 


duF 

~dT 


n 


(4.48) 
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Fluid  mesh/Structure  Interface  conditions 


X  =  US  (4.49) 


dx  _  dus 
dt  dt 


(4.50) 


4.6.2  Numerical  Aeroelastic  Simulations 

The  altitude  was  fixed  at  3,000  m,  the  angle  of  attack  was  set  to  0  degree,  and  the  first  of  a  series 
of  aeroelastic  simulations  was  performed  for  a  sequence  of  stabilized  flight  conditions  at  the  following 
Mach  numbers:  0.8,  0.85,  0.875,  0.9,  0.925,  0.95,  1.0,  1.1,  1.2,  1.3,  and  1.4.  For  each  Mach  number, 
the  numerical  simulation  starts  from  a  steady-state  flow  and  the  following  initial  conditions  for  the 

typical  wing  section:  h°  =  0.01  m/s,  and  6°  -  0.2  rad/s.  Fig.  4.5  and  Fig.  4.6  present  the  predicted 
transient  responses  of  the  structure  at  Mx  =  0.8 ,  and  Mk  =1.0,  respectively.  These  figures  show  that 
at  M x  =  0.8  ,  the  aeroelastic  vibrations  are  rapidly  damped  out.  They  also  show  that  at  Mx  =  0.8  , 
both  modes  of  the  typical  wing  section  contribute  initially  to  the  bending  degree  of  freedom  ( 0 ).  It 
can  be  concluded  that  it  is  the  bending  mode  of  the  typical  wing  section  that  is  rapidly  damped  out  at 
Mx  =  0.8  .  The  history  of  6(t)  graphically  depicted  in  Fig.  4.6  reveals  that  flutter  or  limit  cycle 
oscillations  are  initiated  at  Mx  =  1.0  . 
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Figure  4.5  Transient  aeroelastic  response  for  stabilized  flight  conditions  at  Mach  =  0.80 


Figure  4.6  Transient  aeroelastic  response  for  stabilized  flight  conditions  at  Mach  =  1.00 


Next  the  aeroelastic  response  during  accelerated  flight  of  the  typical  wing  section  is  simulated  at 
the  rate  of  0.05  Mach/second.  The  reader  is  reminded  that  such  acceleration  is  even  higher  than  what 

an  F-16  can  achieve  in  a  level  flight.  This  simulation  was  initiated  at  =  0.75  and  excite  the 

structure  with  the  same  initial  conditions  as  previously:  h°  =  0.01  m/s,  and  0°  =  0.2  rad/s.  It  is 
reported  in  Fig.  4.7  the  predicted  response  of  the  structure.  The  reader  can  observe  that  both  modes  of 
the  typical  wing  section  contribute  to  the  early  response  of  the  h  degree  of  freedom.  The  significant 
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decrease  of  the  mean  value  of  h  that  is  noted  between  Mach  0.85  and  Mach  0.95,  suggests  an 
important  drop  in  the  lift  in  this  Mach  region,  which,  given  that  the  F-16  airfoil  is  unsymmetric.  This  is 
indicative  of  the  appearance  of  a  shock  and  reaching  the  Mach  divergence  speed.  Most  importantly, 
Fig.  4.7  shows  that  between  Mach  0.85  and  Mach  0.95,  the  vibrations  of  h  and  6  become  too  small 
to  allow  a  parametric  identification.  Hence,  it  can  be  concluded  that  a  continuous  parametric 
identification  of  an  aircraft  that  accelerates  across  the  subsonic,  transonic,  and  supersonic  regimes 
requires  systematic  re-excitations. 

In  order  to  illustrate  the  effect  of  re-excitation,  a  second  accelerated  flight  was  simulated  where 
the  initial  speed  corresponds  to  Mach  number  0.84,  and  the  initial  excitation  is  effected  by  the  same 
initial  conditions  as  previously.  The  computed  aeroelastic  response  is  reported  in  Fig.  4.8.  This 
response  is  characterized  by  larger  amplitudes  of  vibration  between  Mach  0.85  and  Mach  0.95,  and 
confirms  that  flutter  or  limit  cycle  oscillations  initiate  around  Mach  1.0. 
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Figure  4.7  Transient  response  during  an  accelerated  flight:  initial  speed  corresponds  to  Mach  0.75  and  is 
increased  by  0.05  Mach/second 
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Figure  4.8  presents  the  simulated  aeroelastic  response  of  the  accelerated  typical  wing  section  after  re¬ 
excitation  at  Mach  1.04. 


Figure  4.9  Transient  response  during  an  accelerated  flight:  initial  speed  corresponds  to  Mach  1.04 
and  is  increased  by  0.05  Mach/second 


4.63  Validation  of  the  2  Dimensional  Simulation 


As  stated  in  Section  4.0,  the  ERA  assumes  a  constant  sampling  rate.  Since  the  numerical 
simulations  are  not  performed  with  a  constant  time-step  At,  the  simulation  results  are  post-processed 
by  using  a  quadratic  interpolation  scheme  in  order  to  generate  signals  with  a  constant  sampling  rate. 
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The  typical  wing  section  described  in  this  chapter  has  2  degrees  of  freedom.  However,  the 
aeroelastic  typical  wing  section  has  more  than  2  degrees  of  freedom,  because  of  the  surrounding  fluid. 
For  this  reason,  in  all  cases  discussed  in  this  section,  the  number  of  states  for  the  synthesized  system  is 

set  to  nx  =  10 ,  as  if  the  system  contained  5  degrees  of  freedom. 

ERA  requires  two  cycles  of  the  lowest  mode  contributing  to  the  signal,  and  about  500  sampling 
points  per  cycle  in  order  to  accurately  identify  the  system.  Hence,  in  this  work  the  parameters  of  the 
ERA  are  set  as  follows 


1000  <V<  6000 
500<<7<1000 
500  <d  <5000 


<At  <3 


ZL 

N 


nx=  1° 


Where  A tt  is  the  sampling  rate,  and  T  is  the  period  of  the  lowest  frequency  mode.  It  was  verified 

a  posteriori  that  the  ERA  configured  with  the  above  parameters  produced  excellent  results  for  all 
applications  discussed  in  this  work.  Nevertheless,  it  was  noted  that  higher  values  of  the  number  of 
samples  N  and  higher  values  of  q  improve  the  accuracy  of  the  identification,  but  increase  its  cost. 

As  explained  in  Section  4.0,  in  order  to  extract  the  aeroelastic  parameters  of  the  typical  wing 
section  from  the  signals  generated  by  the  accelerated  flight  simulations,  ERA  with  windowing  were 
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employed.  The  size  of  each  time-window  is  NMt .  Given  that  the  frequency  of  the  first  wet  mode  of 

the  system  can  be  expected  to  be  close  to  that  of  the  first  dry  mode  of  the  system  -  that  is  4.7  Hz  -it 
follows  that  the  size  of  each  time-window  varies  between  0.4  second  and  0.6  second.  Hence,  for  an 
acceleration  of  0.05  Mach/second,  the  maximum  variation  of  the  Mach  number  within  a  time-window 
is  about  0.03  Mach.  The  significance  of  this  variation  depends  on  how  fast  the  damping  coefficient  of 
the  structure  varies  with  the  Mach  number,  which  depends  on  the  flow  regime.  Extensive  experiments 
have  revealed  that  when  using  the  ERA  with  windowing,  the  identification  results  are  to  some  extent 

insensitive  to  small  variations  in  the  size  of  the  time-window  and/or  the  values  of  flx ,  q ,  and  d  . 

The  frequencies  and  damping  ratios  identified  by  the  FastERA  using  stabilized  and  accelerated 
flight  simulation  data  are  reported  in  Fig.  4.10  and  Fig.  4.1 1.  The  following  observations  are 
noteworthy 

•  The  typical  wing  section  exhibits  flutter  for  the  second  (torsional)  mode  in  the  transonic  regime, 
for  0.89  <  M ^  <1.1. 

•  A  sharp  decrease  of  the  damping  ratio  S2  of  the  torsional  mode  occurs  around  Mach  0.9,  and  a 
slow  increase  of  that  damping  ratio  occurs  above  Mach  0.92.  The  same  trend  was  also  observed  for 
the  damping  ratio  s  of  the  first  mode;  however,  Sx  remains  positive. 

•The  frequency  of  the  first  mode  appears  to  be  almost  independent  of  the  Mach  number.  On  the 
other  hand,  the  frequency  for  the  torsional  mode,  which  is  responsible  here  for  flutter,  increases 
slightly  with  the  Mach  number. 

•  Occasionally,  the  frequencies  and  damping  ratios  identified  for  the  first  mode  in  simulated 
accelerated  flight  are  reported  to  be  different  from  those  identified  using  stabilized  flight  simulated 
data.  This  is  due  to  the  fact  that  in  the  accelerated  flight  simulation,  because  the  system  is  excited  only 
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Mode  1  (bending) 


Mach 


Figure  4.10:  Identified  first  bending  (Mode  1) 
frequencies  and  damping  ratios  for  stabilized  and 
accelerated  simulation  data  from  the  typical  wing  section 
of  an  F- 16  aircraft. 
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Mode  2  (torsional) 


o  from  accelerated  flight  (0.05  Mach/second) 
t  from  stabilized  flight  conditions 

Figure  4.1 1:  Identified  first  torsion  (Mode  2)  frequencies 
and  damping  ratios  for  stabilized  and  accelerated 
simulation  data  from  the  typical  wing  section  of  an  F-16 
aircraft. 


So  it  is  reasonable  to  conclude  that  for  the  simulated  data  the  acceleration  does  not  significantly 
affect  the  aeroelastic  parameters.  It  can  also  be  concluded  that  using  the  windowing  approach  the 
aeroelastic  parameter  could  be  identified  using  the  accelerating  data. 

Now  the  simulated  data  and  the  flight  test  data  provided  by  the  Air  Force  Flight  Test  Center  are 
compared.  The  simulation  results  and  test  data  are  for  the  same  stabilized  flight  conditions,  but  not 
exactly  the  same  F-16  configuration.  The  typical  wing  section  designed  in  this  paper  is  for  a  clean¬ 
wing  configuration  of  the  F-16.  The  flight  test  data  provided  by  the  Air  Force  Flight  Test  Center  are 
for  a  configuration  of  the  F-16  that  includes  pylons  and  missiles.  Nevertheless,  Fig.  4.1 1  shows  that 
the  predicted  wet  frequencies  are  in  good  agreement  with  those  measured  in  flight  test.  However,  the 
predicted  damping  ratios  do  not  agree  well  with  the  experimental  data,  even  though  the  trend  of  their 
variation  with  the  Mach  number  is  similar  to  that  of  the  flight  test  data.  It  can  reasonably  be  argued 
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that  a  major  cause  of  this  discrepancy  is  the  typical  wing  section  model,  which  is  supposed  to  be 
realistic  only  for  fairly  homogeneous  wings  with  high  aspect  ratios  and  angles  of  sweep. 

It  has  been  shown  that  ERA  can  be  used  with  the  windowing  technique  developed  to  identify 
modal  parameter  of  an  accelerating  aeroelastic  system.  It  has  also  been  demonstrated  by  way  of 
simulation,  not  only  that  the  modal  parameters  can  be  identified  from  the  data  of  an  accelerating 
system,  but  the  acceleration  does  not  significantly  change  the  modal  parameters  of  the  system  at  any 
particular  Mach  number.  The  results  of  this  validation  were  so  compelling  that  the  Air  Force  Flight 
Test  Center  agreed  to  conduct  a  test  program  during  which  they  flew  accelerated  and  stabilized  test 
points.  The  details  of  this  test  program  are  the  subject  of  next  chapter,  and  the  results  will  be  used  in 
chapter  6  to  accomplish  further  validation  of  the  accelerated  flutter  test. 


Mode  2  (torsional) 


V  flight  test  data 

Figure  4.12:  Validation  of  ERA  using  windowing  and 
Accelerated  Flutter  Testing  using  pre-1998  data  and 
two  degree  of  freedom  F-16  model 
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5.0  Flight  Test  of  an  F-16  Configuration 

Encouraged  by  the  results  of  the  two  degree  of  freedom  F-16  typical  wing  section  model  the  Air 
Force  Flight  Test  Center  wanted  to  see  results  for  full  aircraft  simulations,  and  supported  the  effort 
with  flight  test  data  flown  specifically  for  this  project  by  the  United  States  Air  Force  Test  Pilots  School 
(TPS).  In  May  of  2000  the  TPS  flew  three  sorties  in  order  to  collect  flight  test  data  for  aeroelastic  and 
Limit  Cycle  Oscillation  conditions.  It  is  interesting  to  note  that  five  sorties  were  planed,  but  all  the 
data  were  collected  in  two  sorties  (one  sortie  was  unproductive  due  to  instrumentation  problems).  One 
explanation  for  this  unexpected  efficiency  is  the  points  were  being  flown  using  accelerations  rather 
than  stabilized  test  points.  The  testing  was  done  as  a  staff  project  of  the  TPS.  The  TPS  was  designated 
the  Responsible  Test  Organization  (RTO)  and  the  416  FLTS  (F-16  Combine  Test  Force)  was 
designated  the  Participating  Test  Organization  (PTO).  This  project  was  named  “Have  Zip”. 


5.1  Test  Item  Description 

The  F-16C  aircraft,  serial  number  S/N  87-0352  was  the  test  aircraft  used  for  these  tests.  The 
aircraft  was  equipped  with  a  flutter  excitation  system  (FES).  This  system  has  the  capability  to  drive 
the  flaperon  servo  actuators  symmetrically  or  antisymmetrically.  The  FES  can  provide  random  burst 
inputs  or  sinusoidal  burst  or  sweep  inputs  at  flaperon  amplitudes  up  to  ±1  degree. 
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5.2  Overall  Test  Objective 

The  test  objective  for  this  flight  testing  was  to  collect  flight  test  data  for  the  F- 16  in  a  clean  wing 
and  a  LCO  configuration  to  support  the  validation  of  the  accelerated  flutter  flight  test  technique. 


5.3  Test  Resources 


5.3.1  Test  Aircraft  Loadings 

Testing  was  performed  using  2  loadings  as  shown  in  table  5.1.  The  aircraft  was  flown  in  both  the 
clean  wing  (no  external  stores)  and  in  a  store  loading  that  provided  moderate  LCO  condition.  Mass 
properties  information  in  the  vertical  and  lateral  directions  were  required  for  all  stores  used  on  this 
program.  The  following  stores  were  required: 

1)  Four  AIM-120  missiles  (AMRAAM) 

2)  Six  LAU-129  launchers 

3)  One  300  gallon  fuel  tank  (Optional) 
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Table  5.1 

Have  Zip  Test  Loadings 


a 

o 

□ 


AIM- 120 


300  Gallon  Tank  on  16S951  Universal  Centerline  Pylon 


LAU-129  Pylon 


5.4  Test  Execution 


The  purpose  of  this  test  program  was  to  acquire  actual  aircraft  accelerometer  and  strain  gage 


response  data.  The  data  are  being  used  to  validate  the  accelerated  flutter  testing  approach. 
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5.4.1  Test  Objectives 

Test  Objective  1 

Collect  baseline  dynamic  response  data  from  an  F-16  Block  40  aircraft  (without  stores) 
during  accelerations  to  supersonic  conditions  and  during  elevated  g  conditions  in  which  the 
aircraft  flexible  modes  are  being  excited  regularly. 

Measure  of  Performance  (MOP)l 

All  clean  loading  test  points  in  Appendix  2  have  been  performed  and  are  flown  within 
specified  tolerances. 

Success  Criteria 

MOP  1  was  satisfied  when  all  critical  data  were  acquired  for  each  test  maneuver. 

Test  Objective  2 

Collect  dynamic  response  data  from  an  F-16  Block  40  aircraft  with  a  store  loading 
yielding  LCO  during  accelerations  to  supersonic  conditions  and  during  elevated  g  conditions 
in  which  the  aircraft  flexible  modes  were  being  excited  regularly. 

MOP  2 

All  store  loading  test  points  in  Appendix  2  were  performed  and  were  flown  within 
specified  tolerances. 


Success  Criteria 


72 


M0P2  was  satisfied  when  all  critical  data  were  acquired  for  each  test  maneuver. 


5.4.2  Data  Requirements 

The  data  for  all  test  points  were  acquired  using  Interactive  Analysis  and  Display  System  (IADS) 
[27].  The  data  were  then  converted  to  engineering  units  in  an  ASCII  format  using  IADS. 


5.5  Test  Procedures 

Flight  testing  consisted  of  two  flight  test  loadings.  One  clean  wing  and  one  store  loading.  The 
clean  wing  loading  was  flown  first  on  12  May  2000  (flight  705).  The  store  loading  was  flown  on  15 
May  2000  (flight  706).  After  flight  705  it  was  discovered  that  the  instrumentation  was  not  functioning 
properly  so  the  clean  wing  maneuvers  were  flown  again  on  16  May  2000  (flight  707).  Test  maneuvers 
consisted  of  straight  and  level  accelerations  from  0.6  to  1 .2  Mach  number  and  windup  turns.  All 
maneuvers  began  at  10,000  feet  Mean  Sea  Level  (MSL).  The  flutter  excitation  system  (FES)  was  used 
in  the  antisymmetric  flaperon  random  burst  mode  for  each  test  point  in  order  to  excite  the  aircraft 
structural  modes.  The  general  test  procedures  were  as  follows. 

1 .  Set  up  at  10,000  feet  MSL  and  0.6  Mach  number  and  perform  level  accelerations  to 

1.2  Mach  number  while  activating  minumum  duration  antisymmetric  random  bursts  (at 
about  1  second  intervals)  throughout  the  runs.  The  excitation  amplitudes  may  be  altered 
in  flight  depending  on  in-flight  response  data.  The  first  level  acceleration  was  performed 
starting  at  MIL  power  (maximum  power  available  without  using  afterburner)  and  then 
throttling  to  MAX  power  (Maximum  available  power  using  afterburner)  at  an 
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acceleration  rate  of  approximately  2-3  knots  per  second.  The  second  level  acceleration 
will  be  performed  using  MAX  power  from  the  beginning  to  the  end  of  the  maneuver. 

2.  Set  up  1 0,000  ft  MSL  and  test  Mach  Number.  Perform  a  Wind  Up  Turn  (WUT)  to  a 

target  load  factor  and  initiate  an  antisymmetric  random  burst  while  maintaining  load 
factor.  Perform  2  more  random  bursts  allowing  response  to  decay  (approximately  3 
seconds)  before  initiating  next  burst. 


5.6  Test  Execution 


5.6.1  Instrumentation  Checkout 

Prior  to  initiating  any  test  maneuvers  on  each  flight.  The  pilot  performed  an  instrumentation 
validation  check.  The  instrumentation  check  was  performed  at  0.8  Mach  number,  10,000  feet  MSL. 
The  loads  instrumentation  validation  check  consisted  of  the  following  maneuvers: 

1 )  3 .5g  windup  turn  (WUT) 

2)  lg  180  degree  roll  in  both  directions 

3)  Wings  level  sideslip  (WLSS)  in  both  directions  with  slow  release. 


5.6.2  Test  Approach 

A  crew  consisting  of  a  pilot  test  conductor,  test  engineer  were  present  for  each  flight.  All  test 
points  were  initiated  at  10,000  ft  MSL.  Test  points  consisted  of  level  accelerations  and  WUT’s  were 
perfomed  at  0.6,  0.75,  0.9,  0.95,  1.05,  1.2  Mach  number. 


74 


5.6.3  Test  Procedures 

Prior  to  setting  up  at  the  test  condition,  the  pilot  would  verify  that  the  appropriate  FES  selections 
were  current.  For  the  level  acceleration  test  points,  and  for  the  WUT’s  the  pilot  stabilized  at  the 
beginning  Mach  number  for  10  seconds.  The  test  engineer  provided  clearance  to  initiate  the  maneuver 
through  the  test  conductor.  The  pilot  then  performed  the  maneuver. 

For  the  clean  configuration,  there  was  no  requirement  in  the  order  in  which  the  test  points  were 
accomplished.  The  desired  test  sequence  was  as  follows: 

1)  Set  up  at  the  test  Mach  number  in  trimmed,  straight  and  level  flight. 

2)  With  the  throttle  in  military  rated  power  (maximum  power  available  without  using 
the  after  burner),  apply  throttle  as  needed  to  begin  a  level  acceleration  to  1 .2  Mach 
number  exercising  the  FES  with  1-2  seconds  between  bursts. 

3)  Reduce  speed  to  0.6  Mach  number  and  stabilize  at  10,000  ft  MSL,  0.6  Mach 
number. 

4)  With  the  throttle  in  Max  power,  begin  a  level  acceleration  to  1 .2  Mach  number 
exercising  the  FES  with  1-2  seconds  between  bursts. 

5)  Reduce  speed  to  0.6  Mach  number  and  stabilize  at  10,000  feet  MSL,  0.6  Mach 
Number 

6)  Initiate  a  left  WUT  and  hold  at  3g. 

7)  Run  a  2-second  antisymmetric  random  burst  three  waiting  3  seconds  between  bursts 
to  allow  the  aircraft  response  to  dampen. 

8)  After  completion  of  the  random  excitation,  pause  for  3  seconds  and  increase  load 
factor  to  5g. 

9)  Repeat  step  7. 
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10)  Pause  for  3  seconds  after  completion  of  random  excitation,  then  return  aircraft  to  1-g 
level  flight  at  the  target  Mach  Number  and  altitude. 

1 1)  Repeat  steps  6-10  for  each  remaining  target  Mach  number  building  up  in  Mach 
Number. 

These  steps  were  the  same  as  were  used  for  the  store  loading  except  that  the  test  point  sequence 
must  follow  the  buildup  sequence  outlined. 


5.6.4  Test  Maneuvers 

5.6.4. 1  Level  Acceleration 

The  aircraft  was  trimmed  to  a  1-g  flight  condition.  A  level  acceleration  was  performed 
to  the  end  Mach  number  initiating  the  acceleration  at  MIL  or  MAX  thrust  as  determined  by 
the  test  point.  The  maneuver  was  ended  at  MAX  thrust  setting.  During  the  acceleration,  2- 
second  antisymmetric  random  bursts  were  performed  at  the  specified  burst  amplitude  with  1-2 
seconds  between  bursts. 

5. 6.4. 2  Wind  Up  Turn 

The  aircraft  was  trimmed  to  a  1-g  flight  condition.  A  stabilized  g  turn  was  performed  by 
banking  the  aircraft  into  a  turn  and  applying  longitudinal  stick  force  to  increase  load  factor  at 
a  rate  of  less  than  2g  per  second  to  the  target  load  factor.  The  throttle  could  be  modulated  as 
required  to  maintain  test  conditions.  The  pilot  then  initiated  a  single  2-second  random  burst  at 
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the  specified  burst  amplitude  and  allowed  the  aircraft  response  to  dampen  for  3  seconds  while 
load  factor  was  mantained.  The  pilot  repeated  the  bursts  until  3  bursts  had  been  performed. 


5. 6. 4. 3  Tolerances 


The  following  tolerances  were  adhered  to  for  all  test  maneuvers: 


Mach  number 

±0.02 

Altitude 

+1000  ft 

Nz 

±0.2g 

5.6.5  Post  Test  Briefing 

After  each  mission  a  post  flight  briefing  was  conducted  to  discuss  aircraft  status,  instrumentation 
status,  review  the  test  points  completed,  and  plan  the  next  flight. 


5.7  Flight  Test  Results 

The  frequency  and  damping  plots  for  the  clean  wing  data  that  were  collected  during  this  test 
program  are  presented  in  Fig  5.1  and  5.2.  Figure  5.1  shows  the  data  for  the  accelerated  flight  test 
points.  The  curve  that  is  designated  as  “normal”  corresponds  to  the  accelerated  test  point  that  was 
flown  using  military  rated  power  (i.e.  maximum  available  power  without  using  afterburner).  The 
curve  that  is  designated  as  “maximum”  corresponds  to  the  accelerated  test  point  that  was  flown  using 
afterburner.  Two  effects  that  contribute  to  the  discrepancy  between  the  two  test  points  are  the 
vibration  environment  caused  by  the  afterburner,  and  the  fuel  configuration.  The  afterburner 
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significantly  increases  the  ambient  vibration  environment,  which  raises  the  noise  floor  for  the  data,  and 
reduced  the  signal-to-noise  ratio.  The  erratic  damping  measurements  are  likely  due  to  this  effect. 
Because  the  test  points  were  performed  at  different  times  during  the  test  flight  there  would  be  a 
different  fuel  load  in  each  case.  A  change  in  the  fuel  load  can  change  the  frequency  and  damping 
characteristics  of  the  system  under  test.  The  change  in  frequency  between  the  two  test  points  is 
possibly  the  result  of  different  fuel  loads. 


Figure  5.1:  Frequency  and  damping  plots  for  flutter  flight  tests  utilizing  stabilized,  accelerated 
normal  thrust,  and  accelerated  maximum  thrust  conditions 


Figure  5.2  presents  the  elevated  g  maneuvers  contrasted  against  the  normal  level  acceleration. 

One  interesting  trend  is  where  the  mach  dip  occurs  in  each  case.  As  the  g  factor  increases,  the  angle  of 
attack  also  has  to  increase.  When  the  angle  of  attack  increases  the  flow  acceleration  over  the  wing,  the 
unsteady  shock  wave  that  causes  the  Mach  dip  will  occur  at  a  lower  Mach  number  consistent  with  the 
flight  test  data  that  was  collected. 


Torsional  Frequency  (Hz) 
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Figure  5.2:  frequency  and  damping  data  for  lg,  3.5g,  and  5g  flight  test  data 


This  flight  test  data  will  be  used  in  chapter  6  in  order  to  validate  the  accelerated  flutter  flight  test 
approach. 
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6*0  3-D  Simulation  of  an  F-16  Configuration 

In  this  chapter  the  concept  of  accelerated  flutter  flight  test,  and  that  aeroelastic  parameters  are 
insensitive  to  acceleration  will  he  further  validated  using  3-D  simulations.  Recall  that  in  chapter  4  the 
validation  using  2-D  simulations  yielded  encouraging  results,  but  still  left  room  for  improvement.  The 
flight  test  data,  for  an  F-16  Block  50  aircraft  (shown  in  Fig  6.1)  with  no  external  stores  (clean  wing) 
that  are  described  in  chapter  5,  and  full  3-D  simulation  data  of  this  F-16  configuration  will  be 
compared  to  perform  this  validation.  An  additional  benefit  to  validating  the  accelerated  flutter  flight 
test  approach  using  simulations  is  that  the  simulation  technology  will  also  be  validated.  The  F-16  was 
selected  for  these  full  aircraft  simulations  because  models  were  obtainable.  Lockheed-Martin,  Fort 
Worth  with  the  approval  of  the  F-16  Systems  Program  Office  released  two  structures  models.  The  Air 
Force  Research  Lab,  Wright  Patterson  AFB  released  the  surface  geometry  that  was  used  to  create  the 
fluid  grid.  These  models  required  many  refinements  in  order  to  bring  them  up  to  the  resolution  needed 
for  the  simulations.  Before  the  validation  data  are  presented,  an  overview  of  the  structural,  and 
aerodynamic  model  refinements,  and  the  simulation  technology  will  be  presented. 
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Figure  6.1:  Air  Force  F-16  Block  50  military  fighter  aircraft 


6.1  Structural  Model  of  the  F-16 


The  models  that  were  obtained  from  Lockheed-Martin  were  a  structural  dynamic  model  (shown  in 
figure  6.2),  and  a  higher  fidelity  static  loads  model  (shown  in  figure  6.3). 


Figure  6.2:  F-16  Block  50  structural  dynamic 
model  Provided  by  Lockheed-Martin. 


Figure  6.3:  F-16  Block  50  static  loads 
model  provided  by  Lockheed-Martin. 
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Due  to  the  computational  resources  that  were  available  when  these  models  were  created,  the 
resolution  of  the  structural  dynamics  only  has  about  1 000  degrees  of  freedom.  Because  this  model  is 
so  coarse  it  was  necessary  to  merge  the  mass  data  with  the  static  loads  model.  In  addition  to  the 
addition  of  mass  to  the  static  loads  model  there  were  a  number  of  refinements  required  so  it  would 
match  the  aerodynamic  model.  These  refinements  were  needed  so  that  the  structural  model  would 
match  the  aerodynamic  model.  It  is  necessary  for  the  structural  model  to  match  the  aerodynamic 
model  so  that  the  software  can  determine  where  on  the  structure  to  apply  the  aerodynamic  (pressure) 
loads.  The  following  figures  show  in  detail  the  refinements  that  were  made  in  order  to  achieve  this 
matching.  Figure  6.4  shows  the  nosecone  and  part  of  the  canopy  missing  in  the  original  static  loads 
model.  These  were  not  modeled  originally  because  they  are  not  structural  members.  However  they 
were  needed  for  the  aeroelastic  simulations.  Figure  6.5  shows  the  final  composite  structural  dynamics 
model  with  phantom  elements  representing  the  nosecone  and  canopy.  Figure  6.6  shows  the  horizontal 
stabilizer  of  the  original  static  loads  model  and  the  overlaid  surface  geometry.  Notice  the  horizontal 
stabilizers  look  like  rectangular  blocks,  and  they  don’t  match  the  anhedral  angle  of  the  actual 
geometry.  Originally  there  was  no  need  to  match  the  geometry  exactly,  but  for  the  aeroelastic 
simulations  it  will  be  necessary  so  that  the  aerodynamic  pressure  loads  could  be  properly  applied  to  the 
structural  surface.  Figure  6.7  shows  the  gap  between  the  control  surfaces  and  the  wing.  This  gap  was 
bridged  using  phantom  elements.  Figure  6.9  show  the  final  composite  structural  dynamics  model 
without  the  gap.  This  was  an  extensive  effort  that  took  approximately  9  man  months. 
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Figure  6.4:  Nosecone  and  att  canopy  missing  on  Figure  6.5:  Phantom  elements  represent 

Original  static  loads  model  nosecone  on  final  model 


Figure  6.6:  Horizontal  stabilizers  do  not 
Match  the  surface  geometry  of  the  aircraft 
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Figure  6.7:  Gaps  between  wing  and  control  Figure  6.8:  Gaps  filled  with  phantom 

surfaces  on  original  static  loads  model  elements  on  final  structural  dynamic  model 


Figure  6.9:  Final  structural  dynamic  model 


6.2  Euler  Fluid  Grid  of  the  F-16 


The  aerodynamic  portion  of  the  simulations  was  computed  using  the  Euler  equations  for  invisid 


fluid  dynamics.  In  order  to  accomplish  these  calculations  it  was  first  necessary  to  develop  an  Euler 
fluid  grid.  The  first  step  toward  accomplishing  this  was  to  obtain  the  Computer  Aided  Design  (CAD) 
surface  geometry  of  an  F-16  Block  50  aircraft.  The  surface  geometry  shown  in  figure  6.10  was 
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provide  by  the  Air  Force  Research  Laboratory,  Wright  Patterson,  with  the  help  of  contacts  at  the  Air 
Force  Flight  Test  Center  at  Edwards  Air  Force  Base  in  California. 


Figure  6.10:  Computer  Aided  Design  (CAD)  Surface  Geometry 
of  an  F-16  Block  50  Fighter  Aircraft 


Once  the  surface  geometry  was  obtained  a  surface  grid  was  developed  as  an  intermediate  step 
toward  the  full  volume  grid.  Figure  6.1 1  shows  the  surface  grid,  and  figure  6.12  shows  the  full  volume 
grid,  which  has  403,919  grid  points.  This  is  the  volume  grid  that  was  used  for  all  of  the  simulations 


results  presented  in  this  document. 
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Figure  6.1 1:  Surface  Grid  for  the  F-16  Block  50  Figure  6.12:  Volume  Grid  for  the  F-16  Block  50 
Fighter  Aircraft  Fighter  Aircraft 


6.3  Aeroelastic  Numerical  Solutions  for  the  3  Dimensional  Model 

Simulations  of  the  full  three-dimensional  F-16  were  conducted  for  Mach  numbers  of  0.7,  0.8,  0.9, 
1.0,  1.1,  1.2,  1.3,  and  1.4  at  pressure  and  density  values  equivalent  to  3000  M  (10,000  ft).  In  order  to 
validate  the  accelerated  flutter  flight  testing  approach  accelerated  test  point  at  the  same  conditions  and 
at  an  acceleration  rate  on  0.05  Mach/sec.  Simulations  were  also  conducted  at  elevated  3.5  and  5  g’s  in 
the  normal  direction  in  order  investigate  the  effects  of  normal  acceleration.  The  frequency  and 
damping  results  of  the  simulation  for  the  one  g  maneuvers  was  compared  to  the  frequency  and 
damping  results  for  the  flight  test  data  shown  in  figure  6.13. 
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“  3'  3D  Simulation  of  Stabilized  Flight  (Clean  Wing) 

[— !  “*■  3D  Simulation  of  Accelerated  Flight  (Gean  Wing) 

^  Accelerated  Flight  Test  Data  (Gean  Wing,  Normal,  2000) 
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Figure  6.13  shows  the  correlation  between  the  flight  test  data  and  the  simulation  data. 


The  frequency  and  damping  results  for  the  3.5g  and  5g  flight  test  data  were  compared  to 
frequency  and  damping  results  for  the  simulation  data.  Because  the  3.5g  and  5g  flight  test  data  were 
collected  using  the  “wind  up  turn”  technique  (which  is  literally  flying  around  in  circles)  the  aircraft 
was  flying  into  its  own  turbulence  wake.  As  a  result  the  data  are  somewhat  noisy,  and  difficult  to 
identify  the  frequency  and  damping.  Consequently,  one  cannot  expect  the  high  g  data  to  correlate  as 
well  with  simulation  data  as  was  seen  in  the  one  g  case.  Although  interpretation  of  the  data  are 
complicated  by  the  turbulence  encountered  in  the  flight  testing,  the  trend  of  the  data  should  be 
conserved.  The  data  shown  in  figure  6.14  compares  the  frequency  and  damping  values  for  the  3.5g 
flight  test  data  with  frequency  and  damping  values  for  the  3.5g  simulation  data.  As  one  can  observe, 
the  frequency  data  correlates  well,  and  the  trend  of  the  damping  values  is  captured. 
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Figure  6.14:  frequency  and  damping  correlation  between  flight  test  data  and  simulation  data  for  a  3.5g 
maneuver. 


In  figure  6.15  the  frequency  and  damping  plots  for  the  5g  flight  test  data  are  compared  to  the 
frequency  and  damping  plots  for  the  5g  simulation  results.  Again  the  frequency  plots  agree  very  well, 
while  the  damping  plots  capture  the  trend. 
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Figure  6.15:  frequency  and  damping  correlation  between  flight  test  data  and  simulation  data  for  5g. 


Therefore,  in  every  case,  the  aeroelastic  parameters  are  insensitive  to  acceleration  whether  it  is  in 
the  longitudinal  or  normal  direction.  It  has  also  been  shown  that  using  the  windowing  technique  and 
FastERA  it  is  possible  to  identify  the  aeroelastic  parameters  from  accelerating  data,  whether  it  comes 
from  simulation  or  from  flight  test.  It  can  also  be  concluded  that  for  an  F-16  fighter  configuration 
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without  external  stores  the  simulations  do  a  very  good  job  of  capturing  the  aeroelastic  behavior  of  the 
actual  aircraft  even  though  physical  phenomenon  like  viscosity  and  turbulence  were  not  accounted  for 
in  the  simulation.  These  phenomena  would  become  more  influential  (especially  viscosity)  when 
simulating  a  configuration  with  external  stores  on  the  wings  of  the  aircraft,  or  at  higher  angles  of 


attack. 
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7.0  Real  Time  Predictive  Analysis 

In  this  chapter,  the  third  issue  is  examined,  namely  “the  need  to  anticipate  the  vibration  response 
at  airspeeds  beyond  which  the  aircraft  has  accelerated”.  In  order  to  accomplish  this  an  analytical  study 
of  an  aeroelastic  system  will  be  performed.  It  would  be  untenable  to  consider  an  analytical  study  on 
the  full  aeroelastic  system  of  the  F-16  that  has  been  presented  in  chapter  6.  Therefore  a  simplified 
aeroelastic  system  will  be  considered  for  this  study.  These  simplified  models  will  only  give  an 
indicator  of  the  stability  of  the  system.  They  lack  the  fidelity  to  provide  reliable  predictive  analysis. 
However  what  is  needed  to  support  knock-it-off  decisions  is  an  educated  guess.  As  the  flight  test 
progresses  the  simplified  models  can  be  adjusted  using  feedback  from  the  collected  data.  The  adjusted 
model  will  provide  an  improved  estimate  for  the  stability  of  the  system.  Before  the  aircraft  reaches  the 
critical  speed,  the  model  should  be  refined  enough  to  anticipate  the  danger,  and  the  flight  test  can  be 
halted.  Historical  results  from  flutter  flight  testing  reveals  that  the  flutter  mechanisms  on  full  scale 
aircraft  have  typically  been  dominated  by  one  or  two  modes.  Therefore  only  systems  with  one  or  two 
degrees  of  freedom  will  be  considered  in  this  study.  These  analytical  studies  will  facilitate 
generalizations  that  can  be  made  to  the  behavior  of  these  systems.  The  results  of  this  study  can  be  used 
to  develop  rational  criteria  for  when  testing  should  be  terminated. 
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7.1  The  One  Degree  of  Freedom  System 

The  one  degree  of  freedom  system  considered  here  is  a  typical  two  dimensional  air  foil  that  is  free 
to  move  in  rotation  (equivalently:  pitch  or  torsion)  but  is  restrained  from  moving  in  plunge 
(equivalently:  bending)  as  shown  in  figure  7.L 


Fig  7.1:  two  dimensional  airfoil  with  one  (rotational)  degree 
Of  freedom. 


7.1.1  N  omenclature 

b-Vi  chord 

a  -  position  of  the  center  of  axis  measured  from  the  simi-chord  and  divided  by  b  (see  fig.  7.2) 
U  -  flight  speed 
h  -  bending  displacement 
a  -  rotational  displacement 

Ao  -  density 

S  -  plan  area  =  2b  for  unit  length  wing  section 
ka  -  torsional  bending  stiffness 
(p  -  indicial  function 
Ia  -  rotational  moment  of  inertia 


M  -  Aerodynamic  Moment  acting  on  airfoil 
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Fig  7.2:  a  is  a  nondimensional  measure  of  the  center 
of  elasticity  from  the  center  chord  such  that  -1  indicates 
the  leading  edge  and  1  indicates  the  trailing  edge 


7.1.2  The  One  Degree  of  Freedom  Mathematical  Model 

From  Equation  6.143  on  page  259  of  Principles  of  Aeroelaticity  [2]  the  equation  of  motion  for  a 
two  dimensional  airfoil  allowed  to  move  only  in  rotation  as  shown  in  figure  7.2  is  as  follows. 


1-  dM  • 

dM  ] 

a - 7- a+ 
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da 

CL 

da 

a  =  0 


(7.1) 


Also  from  Principles  of  Aeroelasticity  [2]:  equation  4.171  on  page  120,  the  expression  for  the 
aerodynamic  moment  acting  on  a  two  dimensional  airfoil  allowed  to  rotate  and  plunge  is  as  follows. 
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Where  s  is  the  distance  in  semichords  traveled  since  t  =  0  .  Specifically,  s  =  — ,  and  (p  is  the 

b 

indicial  function  which  is  a  dimensionless  aerodynamic  operator.  This  operator  is  not  consistently 
defined  in  the  literature,  however  2n  will  cause  the  indicial  function  to  asymptotically  approach  unity. 
(p  has  many  forms  depending  on  the  flow  regime,  and  the  geometry  of  the  aircraft.  For  more 
information  on  the  indicial  function  refer  to  page  120  of  reference  [2]. 

It  is  assumed  that  the  bending  degree  of  freedom  will  be  constrained  to  zero  since  only  the 
torsional  degree  of  freedom  is  allowed  to  move.  This  constraint  is  shown  mathematically  in  equation 

(7.3). 


h  =  h  =  h  =  0 


(7.3) 


Combining  equations  (7.2)  and  (7.3)  yields  the  following  equation. 


(7.4) 
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Taking  the  derivative  of  the  aerodynamic  moment  (equation  (7.4))  with  respect  to  rate  of  rotation 
of  the  airfoil  (a)  and  taking  into  account  that  2, 71 ,  p ^  ,b,  and  U  are  all  constant  with  respect  to 

a  yields  equation  (7.5). 
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Similarly  taking  the  derivative  with  respect  to  CC  yields  equation  (7.6). 
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In  summary  equation  (7. 1)  can  be  written  as: 


C,  a+  C2  a+  C3a  =  0 


(7.7) 


Where 
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Note  that  C2  and  C3  are  dependent  on  s  which  is  itself  dependent  on  time.  Therefore  these  are 
time  dependent  coefficients. 


Now,  three  particular  cases  will  be  considered  where  the  indicial  function  will  be  set  to  zero,  a 
constant,  and  an  arbitrary  function  respectively. 
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Case  I,  where  cp  is  zero,  will  be  considered  because  it  simplifies  equation  (7.7)  making  it  a  linear 
second  order  differential  equation  with  constant  coefficients,  Studying  this  simple  case  will,  hopefully 
give  insight  into  the  more  complex  case.  Case  II,  where  (p  is  constant,  is  a  good  approximation  for 
supersonic  flow  with  a  thin  airfoil.  This  case  is  being  investigated  because  it  includes  real  world 
complexity  and  therefore  provides  a  more  credible  tool  for  extrapolating  damping  trends.  Case  III  is 
clearly  the  most  useful  model,  however  it  also  the  most  complex,  and  will  have  to  be  solved 
numerically. 

7.1.2. 1  Case  I:  (p  ~  0 

In  case  I,  where  (p  is  equal  to  zero,  C2  and  C3  can  be  simplified  as  follows. 

C2=7vpxUb\±-a) 

Q=K 

Notice  here  that  C2  and  C3  are  no  longer  dependent  on  time. 

Since  Cx  is  non-zero  equation  (7.7)  can  be  rewritten  as: 

(Z+CzGC+CsGC  =  0  (7.8) 


Where: 
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Cl 


c, 


and 


c3 


Q 

c, 


There  are  three  sub-cases  that  must  be  considered.  Sub-case  I  where 


C2-4Cs  >0 


Sub-case  II  where 


~  2 

C2-4C3  =  0 


Sub-case  III  where 


~  2 

C2-4C3  <0 


For  sub-case  I  the  solution  to  equation  (7.8)  is: 


a  =  B,e*+B2e* 


Where  \  and  /L,  are  defined  as  follows. 
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^1,2  — 
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—  C2  +  'Si  C2— 4C3 


If  the  initial  conditions  are  a(t0)  =  a0  and  a(t0)  =  a  0  it  can  be  shown  that: 


1_ 


/l2a0  -  ao 

f)9  = - 

2  ^2-4 


For  sub-case  II  the  solution  for  equation  (7.8)  is: 


a  =  BxeXt  +  B2teXt  =  (Bx  +  52/) 


Where: 


A,=^=2  =  -^. 


5,  =a0 


B2  -ao-Aa0 


For  sub-case  III  the  solution  for  equation  (7.8)  is: 
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(Bx  cos  at  +  B2  sin  at) 


Where: 


a  = 


Jc2-4C, 

2 


B ,  =«0 


52  —  q:0  +  ^ 

The  term  COS  COt  +  i?2  sin  C Ot )  is  bounded  so  the  stability,  in  this  case,  depends  on  the 


e  2  term.  The  stability  of  this  term  also  depends  on  the  sign  of  Ci . 

At 

Reference  [17]  shows  the  term  e  2  in  the  form  e  ^  where  £  is  the  damping,  and  COn  is  the 
natural  frequency.  This  result  will  come  in  handy  later  in  this  chapter. 

It  is  also  worth  noting  that,  in  the  case  that  (p  is  zero,  \A^—C 3  and  Al  +  A^  =  —  C2  .  In  other 
words,  if  C3  is  positive  then  the  sign  of  A^  will  be  the  same  as  the  sign  of  A^  and  will  be  opposite 
sign  of  Ci .  One  can  see  that  C3  is  positive,  so  A^  will  have  the  same  sign  as  A^  •  The  sign  of  \ 
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and  X2  will  be  opposite  the  sign  of  Ci  which  can  be  determined  by  the  term  —  —  a .  Therefore  one 

can  conclude  that  \  and  ^  are  negative  if  and  only  if  a  <  — .  Therefore  the  airfoil  will  be  stable  if 

3 

and  only  if  the  center  of  elasticity  is  forward  of  the  —  chord. 

So  for  case  1  ( (p  —  0 )  the  solution  and  therefore  the  frequency  and  damping  values  are  directly 
attainable.  Now  consider  case  2  where  (p  =  constant  which  will  give  results  that  are  more 
representative  of  actual  aeroelastic  systems. 


7. 1.2.2  Case  II:  (p  =  constant 

This  case  is  of  limited  practical  use  because  (p  is  only  a  constant  at  supersonic  speeds,  but  it  is  a 
real  physical  case  and  can  give  insight  to  the  behavior  of  aeroelastic  systems  in  general. 
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Notice  here  that  C2  is  dependent  on  time  but  C3  is  not.  This  leads  to  an  equation  of  the  form: 
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a+  (d,  +  d2s )  a+  d2a  =  0 


(7.9) 


Where: 


xp^U 


<*,=■ 


-a 


C, 


(7.10) 


27ipJJ2b2 
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C, 


(7.11) 


and 


J3=-^  (7.12) 

The  obvious  approach  to  getting  a  solution  for  equation  (7.9)  is  to  use  a  series  solution.  This 
approach  was  investigated  and  found  to  not  lead  to  a  closed  form  solution.  The  next  approach  that  was 
tried  was  reformulating  the  problem  as  a  system  of  first  order  equations,  and  finding  an  estimator  on 
a  that  would  bound  the  actual  amplitude. 

Rearranging  equation  (7.9)  yields: 

a  =  -~-(dx  +d2s)a-d3a  (7.13) 


Next  define 
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z,  -  a 


(7.14) 


z2  =  a 


Taking  the  derivative  of  equations  (7.14)  and  combining  with  equation  (7.13)  yields: 


z2  =  a  =  t/jZj 


-(Jj  +t/25)z2 


(7.15) 


Now  define 


Z  = 


Z2 


(7.16) 


Then  rewrite  equations  (7.15)  as: 


0  1 

~d-i  ~(dl  +  d 2) 


(7.17) 


Next  define: 


w  =  (i\ 


+  z207 


(7.18) 


Taking  the  derivative  of  equation  (7.18)  gives: 
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♦  1  _  1  /  •  • 

w  =  — (z2  +  z2)  /2(2z,  z,  +  2z2  z2) 


Z|  Zj  H-  Z2  Z2 

w 


Recalling  equations  (7.15)  it  can  be  seen  that: 


W  = 


zxz2  +  z2(-d3z1  -((ij  +  d2s)z2)  (1  -d3)zxz2  -(dx  +  d2s)z: 


w 


w 


Consider  the  term 


(1  d3)zxz2  (dx  Jrd2s)z2 


It  is  easy  to  see  that: 


(1  d3)zxz2  [dx  Jrd2s}z2  <  |(1  d3)zxz2  {dx-\-d2s}z2 


From  equation  (7.21)  and  the  triangle  inequality  one  can  see  that: 


|(1- d3)zxz2  -{dx  +d2s)z\  <|(l-rf3)||z1|  z2|  +  |(4  +  d2s)|z2 


(hl-hl)2^0 

z\  -2|z,z2|  +  z2  >  0 


z,2  +  z2  >  2  z,z2 


(7.19) 


(7.20) 


(7.21) 


(7.22) 


Now  consider  the  quantity: 
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This  result  and  inequality  (7.22)  to  conclude  leads  to  the  conclusion  that: 


(7.23) 

\{\-d,)\\zx\\z2\  +  \(d,+d2s)\z22< 


H 
2 


■  + 


|(4+4^)| 


w 


Combining  the  inequalities  (7.21)  (7.22)  (7.23)  with  equation  (7.20): 


w<k(s)w 


Where 


k(s)  = 


— y^+|(4  +4>0| 


(7.24) 


Now  apply  Gronwall’s  lemma  [15]. 


s 

J  k(s)dr 

w|  < \wa\es<> 


(7.25) 


Where 


W0=w(s0) 


Recalling  equations  (7.18)  and  (7.24)  equation  (7.25)  becomes: 
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(*T+*J) 


2^K 


< 


(Z12(50)  +  Z22(50))/ 


Recalling  equation  (7.14)  gives: 


jf^7^+l(rfi+rf2s)lVr 

e50 


From  this  it  can  be  deduced  that: 


(a2  +a2) 


Yi 


.  2 


(a0  +ao)7 


a  < 


(#0  +CC  o) 


X 


jflL^l+|(rfl+rf25)|Lr 


{126) 


Equation  (7.26)  is  a  bound  on  the  solution  to  equation  (7.9),  which  represents  the  system  that  is 
under  investigation. 

It  would  be  beneficial  to  have  some  estimate  of  the  damping  of  the  system,  but  the  estimation 
gives  us  a  bound  on  OL  .  In  order  to  get  an  estimate  of  damping  it  is  necessary  to  make  some 
assumptions,  because  damping  is  really  associated  with  a  linear  system,  and  the  system  under 
investigation  is  nonlinear.  However,  an  exponential  bound  has  been  developed  for  the  system: 


es°K  J 


That  resembles  the  exponential  bound  seen  in  case  I,  sub-case  III.  Recalling  that  reference  [17] 


gives  the  solution  in  the  form. 
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e** 

So  if  the  exponent  is  divided  by  —0)J  all  that  remains  is  the  damping  term.  If  it  is  further 

assumed  that  COn  =  0)  where  CO  is  the  damped  frequency  then  the  following  formula  will  give  us  an 
estimate  of  damping. 


l(^+|(rfl+<,2S)|V 

-so'"  ' 

CO 

It  turns  out  that  this  approach  is  extremely  conservative.  That  is  because  the  system  that  is  being 
considered  can  be  divergent,  stable,  or  convergent,  and  the  estimator  implicitly  assumes  the  worst  case, 
that  is  that  the  system  is  divergent. 

A  numerical  study  is  performed  in  section  7.4  that  shows  this  is  indeed  the  case. 


7. 1.2.3  Case  III:  (p  is  arbitrary 

In  the  case  that  (p  is  arbitrary  any  form  of  the  indicial  function  can  be  used.  With  increased 

complexity  an  analytical  solution  becomes  harder  to  come  by.  However,  a  numerical  solution  can 
always  be  reached,  and  this  numerical  solution  can  be  used  to  anticipate  the  damping  at  airspeeds  not 
yet  reached,  and  aid  in  making  the  knock-it-off  decisions  that  sometimes  have  to  be  made  during 
flutter  flight  testing. 
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7. 1 .3  Numerical  Study  of  the  One  Degree  of  Freedom  Error  Analysis 

In  order  to  assess  the  utility  of  the  estimator  (equation  (7.26))  for  the  one-degree  of  freedom 
system,  a  MATLAB  SIMULINK  model  was  built  (shown  in  figure  7.3)  that  includes  the  one-degree  of 
freedom  system,  the  estimator,  and  a  comparison  between  the  two.  The  green  portion  of  the  model  is 
the  one-degree  of  freedom  system  that  will  numerically  solve  equation  (7.9).  The  blue  portion  of  the 
model  is  the  estimator,  which  calculates  the  right  side  of  equation  (7.26).  The  red  portion  of  the  model 
simply  calculates  the  ratio  of  the  outcome  of  the  one-degree  of  freedom  system  and  the  estimator. 
Results  are  presented  for  two  cases  in  figures  7.4  and  7.5.  The  parameters  used  for  each  case  can  be 
found  in  Table  7.1. 

The  first  case  represents  a  set  of  parameters  that  are  not  physically  relevant  and  include  a  very 
small  chord  length  and  very  low  airspeed.  These  parameters  were  selected  such  that  the  estimator 
would  perform  well.  The  data  for  case  I  are  shown  in  figure  7.4.  One  can  see  from  the  figure  that  the 
stability  of  the  system  is  driven  by  the  location  of  the  center  of  elasticity.  If  the  center  of  elasticity  is 
forward  of  the  quarter  chord  the  system  is  unstable.  If  the  center  of  elasticity  is  aft  of  the  quarter  chord 
the  system  is  stable.  This  is  driven  by  d2  term  described  in  equation  (7.1 1).  Specifically  if  the  term 
(!4+a)  is  negative  the  system  will  be  stable.  It  can  also  be  seen  in  the  graph  that  the  estimator  predicts 
the  system  will  be  unstable,  and  the  only  time  the  estimator  is  very  close  to  the  simulation  is  when  a  =  - 
Vi .  One  final  observation  is  that  as  the  airspeed  (U)  increases  the  estimator  becomes  increasingly 
more  conservative. 

The  second  case  represents  a  more  realistic  geometry,  however,  as  one  might  expect  from  the 
trends  observed  in  the  first  case,  the  estimator  is  ridiculously  conservative  when  the  chord  length  and 
the  airspeed  are  more  realistic.  The  second  case  is  illustrated  in  figure  7.5. 

One  positive  thing  that  can  be  said  for  the  estimator  is  that  in  all  cases  it  bounds  the  simulation 


data  perfectly. 
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ConSant3  Producl3 


Display 


1/massl  Product  Stiffness  Displayl 


Figure  7.3:  Matlab  Simulink  model  used  to  simulate  the  one  degree  of  freedom  system  (green), 
estimator  (blue),  and  the  comparison  between  the  two  (red). 
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Table  7.1:  Parameter  values  for  two  cases  presented  for  1  degree  of  freedom  error  analysis  numerical  study. 


Stable  Good  Estimator 

Estimator  for  a  Realistic  System 

Parameter 

Case  1 

Case  2 

U  (mph) 

10  through  100 

100  through  300 

Rho 

0.001 

0.001736 

Ialpha 

0.001 

0.5522 

B 

0.005 

1.0 

S 

0.01 

2.0 

A 

0.0  through -1.0 

-0.49  through  -0.56 

K_alpha 

0.001 

345.1 

Gamma 

1.4 

1.4 

Pressure 

0.001 

1455.7 

aerodynamic  damping  (%) 
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Estimator  performance  for  low  airspeed,  and  small  airfoil 


Figure  7.4:  Comparison  of  damping  results  for  the  estimator  and  simulation  of  the  one  degree  of 
freedom  airfoil  with  a  small  chord  length  and  low  airspeeds. 
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Estimator  performance  for  high  airspeed,  and  large  airfoil 


center  of  elasticity  (a) 


Figure  7.5:  Comparison  of  damping  results  for  the  estimator  and  simulation  of  the  one  degree  of 
freedom  airfoil  with  a  larger  chord  length  and  higher  airspeeds. 


7.1.4  Symbolic  Solution  to  the  One  Degree  of  Freedom  System 


A  solution  to  equation  (7.13)  can  be  found  using  a  symbolic  math  solver.  The  symbolic  math 
solver  in  Matlab  was  used  for  this  purpose.  The  resulting  solution  is  shown  in  equation. 
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Where  W  and  M  are  Whittaker  functions.  Whittaker  functions  are  confluent  hypergeometric 
functions  that  are  defined  below  [16]. 


M =  e  //2zm+//2lFl  ^y^  +  m-k,l  +  2m9z^  (7.28) 

W(k,myz)  =  e  zm+ m-k,\  +  2m,z^  (7.29) 


Where  l  F}  and  U  are  also  confluent  hypergeometric  functions  known  as  Rummer’s  functions. 
Rummer’s  functions  are  defined  below  [16]. 


lF](a;b;z)  =  \  +  ?- 
b 


z  + 


a  + 1)  z 2 
b(b  + 1)  2! 


+ 


(7.30) 


Where  (a)A  and  (b)k  are  Pochhammer  symbols  [16]. 


U(a,b,z) 


r  (a) 


f ez,f-'  (1  +  t)Mdt 


(7.31) 


Where  is  the  gamma  function  [16],  which  extends  the  factorial  function  to  complex  and 


real  numbers. 


Combining  equations  (7.27)  through  (7.31)  results  in  the  following: 


Ill 


-i/4/(^2/+2^) 


a  = 


/^1  + 


C^z"'^  1 


-4-  fe^vri+z)'’  "  '  dt  +  C2e'2z 

!»oJ  V  ' 


/4  vm+^2 


,  01 

1+-Z+- 


a\  a+ 1 


6  6  6+1 


2! 


-+... 


Where: 


*=- 

4 


f 


~  cl 2  +  2  r/3 


V  “2  V 


m=  — 

4 


1 


^(4+</£)2^ 


~  1  ~  :  t  1  d3 

a  =  -  +  m-k  =  1 - - 

2  2d, 


and 


~  ~  3 

b  =  l  +  2m  =  — 

2 

So  given  an  exact  solution  to  the  one  degree  of  freedom  typical  wing  section  it  will  be  possible  to 
project  the  motion  of  the  airfoil  to  any  airspeed  or  Mach  number  as  needed  to  support  knock-it-off 
decisions  in  the  control  room.  A  word  of  caution  is  in  order  here;  this  solution  is  based  on  a  very 
simple  model  and  should  not  be  expected  to  yield  results  that  are  predictive  in  nature.  This  is  only 
intended  to  yield  ballpark  estimates  that  can  be  used  to  support  real-time  decision  making. 


112 


7.2  Knock-it-off  Criteria  &  Real  Time  Displays 

The  final  issue  that  needs  to  be  considered  is  the  knock-it-off  criteria.  Figures  7.6  through  7.8 
depict  a  conceptual  scenario  for  accelerated  flutter  testing.  In  figure  7.6  the  black  predicted  damping 
curve  would  come  from  the  high  fidelity  simulations  that  would  be  conducted  prior  to  the  flight 
testing.  The  red  curve,  labeled  “predicted  damping  adjusted  for  uncertainty”,  is  shown  here  as  simply 
transposition  of  the  original  predicted  curve.  In  principle,  it  could  be  a  more  complex  adjustment  of 
the  original  predicted  damping  curve  that  accounts  for  the  measured  damping  values  that  are  collected 
during  the  flight,  and  the  uncertainties  that  are  present  in  all  such  tests.  The  red  line  at  the  bottom  of 
the  chart  that  is  labeled  “critical  damping”  represents  the  damping  value  that  the  program  will  not 
break  intentionally.  Because  the  parameter  identification  of  data  from  an  accelerating  aircraft  will  add 
an  additional  uncertainty,  a  new  line  labeled  “accelerated  critical  damping”  will  be  used  as  a  “must  not 
break”  line.  The  left  edge  of  the  blue  vertical  band  represents  the  Mach  number  at  the  current  moment 
in  time.  Conceptually  this  band  will  move  across  the  display  from  left  to  right  as  the  aircraft 
accelerates.  The  right  edge  of  this  blue  band  represents  the  expected  Mach  number  that  will  occur 
three  seconds  from  the  current  time.  The  width  of  the  blue  band  will  be  wider  at  higher  acceleration 
rates  and  narrower  at  lower  acceleration  rates.  The  blue  tangent  line  that  is  shown  represents  a 
projection  of  the  damping  trend  that  can  lead  to  a  projected  flutter  speed.  This  could  be  accomplished 
using  the  tangent  line  as  shown,  or  other  technique  such  as  the  flutter-o-meter  concept  developed  at 
NASA  Dryden  [18],  or  the  Zimmerman  method  [19]. 

In  figure  7.7  the  projected  flutter  speed  has  reached  right  hand  side  of  the  band  representing  the 
anticipated  Mach  number  three  seconds  from  the  present  time.  The  band  has  turned  yellow  to 
emphasize  the  problem.  At  this  point  a  real  potential  exists  that  the  aircraft  could  exceed  the  flutter 
speed  before  data  can  be  analyzed  and  control  room  engineers  and  aircraft  crew  can  react  to  avoid  the 


Damping 
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disaster  [20,21].  One  solution  to  this  dilemma  is  to  simply  reduce  the  acceleration  rate  as  shown  in 
figure  7.8.  In  an  extreme  case  the  acceleration  rate  could  be  reduced  to  zero,  in  which  case  the 
accelerated  flutter  test  technique  would  then  degenerate  into  the  conventional  stabilized  flutter  test 
technique. 


Mach 


Figure  7.6:  Proposed  accelerated  flutter  testing  display,  and  methodology. 


Damping  Damping 


Figure  7.8 
or  halted. 
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The  hope  in  conducting  the  research  in  sections  7.1  and  7.2  was  to  develop  a  method  by  which 
damping  could  be  predicted  at  Mach  numbers  higher  than  had  been  flown,  that  would  be  reasonable  (in 
the  ball  park)  and  conservative.  These  predictions  could  then  be  used  to  develop  the  adjusted 
predictive  analysis,  or  perhaps  the  projected  flutter  speed  to  be  used  in  the  knock-it-off  criteria.  For 
the  purpose  of  this  section  we  will  focus  on  the  one  degree  of  freedom  equations  developed  in  section 
7.1  because  the  flutter  mechanism  in  the  F- 16  data  that  we  have  examined  so  far  is  in  fact  a  one  degree 
of  freedom  flutter  mechanism.  In  figure  7.9  the  numerical  solution  of  equation  (7.9),  which  is 
equivalent  to  the  exact  solution  presented  in  equation  (7.27)  is  compared  to  the  two  degree  of  freedom 
model  of  the  F-16  developed  in  chapter  4.  Notice  that  at  the  high  Mach  number  the  system  does  a  very 
good  job  of  predicting  the  damping  of  the  system,  and  recall  that  the  indicial  function  used  to  develop 
equation  (7.9)  was  for  the  supersonic  case.  When  compared  to  the  3  dimensional  flight  test  results 
(7.9)  it  can  be  seen  for  this  case  the  indicial  function  estimator  agrees  well  with  the  flight  test  for  the 
supersonic  condition,  and  is  conservative  for  the  lower  airspeeds. 


2  Dimensional  Comparison  of  Indicial  Function  vs 
Euler  Simulation 


—2  DOF  +  indicia! 
function 

- 2  DOF  +  Euler 


Mach  Number 


Figure  7.9:  Frequency  and  Damping  Comparison  Between  Indicial  Function  and  Euler 
Simulation  for  a  2  Degree  of  Freedom  Typical  Wing  Section  Representation  of  the  F-16. 
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Frequency  Plot:  Indicial  Function  vs 

Flight  Test  vs  3-D  Simulation 
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Figure  7.10:  Frequency  and  Damping  Comparison  Between  Indicial  Function,  3-Dimensional 
Simulation,  and  Flight  Test  Data  for  the  F-16. 
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8.0  Conclusions 

Section  1.6  outlined  several  issues  that  would  need  to  be  addressed  in  order  to  develop  the 
accelerated  flutter  test  technique.  For  convenience  these  issues  are  summarized  below. 

1)  The  first  issue  is  whether  or  not  the  frequencies  and  dampings  are  affected  by 
the  acceleration  of  the  aircraft. 

2)  The  second  issue  is  whether  or  not  parameter  identification  algorithms  can 
extract  frequency  and  damping  values  from  the  time  varying  data. 

3)  The  third  issue  is  the  need  to  anticipate  the  vibration  response  at  airspeeds  (or 
Mach  numbers)  beyond  which  the  aircraft  has  accelerated. 

4)  The  fourth  issue  is  developing  formal  criteria  to  determine  when  the  aircraft 
needs  to  end  the  acceleration  and  terminate  the  test  point. 

The  first  of  these  issues  is  addressed  in  chapter  3.  In  that  chapter  it  was  shown  that  there  are 
inertial  effects  that  mathematically  change  the  basic  aeroelastic  system  that  is  being  investigated. 

These  inertial  effects  change  not  only  the  forcing  function  acting  on  the  aircraft,  but  also  change  the 
stiffness  of  the  torsional  degree  of  freedom  of  the  typical  wing  section.  At  first  glance  this  would 
appear  to  invalidate  the  accelerated  flutter  testing  approach,  however  for  typical  aircraft  geometries 
and  performance  the  change  in  the  torsional  stiffness  (and  therefore  the  flutter  parameters)  will  be  very 
small  (less  than  0.2%  for  the  F-16  case  investigated).  Specifically,  the  position  of  both  the  center  of 
gravity  and  the  center  of  elasticity  are  normally  contained  within  the  wing  of  the  aircraft  and  within 
about  one  quarter  chord  of  each  other.  This  limits  the  moment  arm  that  the  inertial  forces  due  to 
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acceleration  can  act  upon.  Also,  for  level  acceleration,  the  acceleration  of  a  high  performance  aircraft 
is  likely  less  than  0.01  Mach/sec,  which  limits  the  magnitude  of  the  inertial  force  due  to  the 
acceleration  of  the  aircraft.  Because  the  change  in  flutter  parameters  is  negligible  it  is  reasonable  to 
expect  the  accelerated  flutter  test  procedures  to  yield  essentially  the  same  results  as  the  traditional 
stabilized  flight  test  technique,  assuming  the  parameters  can  be  accurately  identified  using  the  data 
from  an  accelerating  aircraft. 

Therefore  it  can  be  concluded  that  the  effect  of  acceleration  on  the  flutter  parameters  is 
negligible  for  normal  geometries  and  acceleration  rates. 

The  second  issue  is  addressed  in  chapter  4.  In  that  chapter  it  was  shown  that  the  wet  frequencies 
and  damping  ratios  of  an  aircraft  are  functions  of  the  free- stream  Mach  number.  Therefore,  during  an 
accelerated  flight,  these  parameters  become  time-dependent.  For  this  reason,  theory  suggests  that 
standard  signal  processing  techniques,  which  are  limited  to  time-invariant  systems,  cannot  be  applied 
to  the  identification  of  an  accelerated  aeroelastic  system.  However,  within  a  time-window  of  0.2 
second,  and  for  a  typical  level  flight  acceleration  of  1  g  or  0.03  Mach/second,  the  Mach  number  varies 
by  0.6%  only.  Hence,  within  a  time-window  of  0.2  second,  the  aeroelastic  parameters  of  an 
accelerating  aircraft  can  be  assumed  to  remain  constant.  Such  a  time-window  corresponds  to  2  cycles 
of  a  mode  at  10  Hz,  a  frequency  that  is  relevant  to  the  first  bending  and  torsion  modes  of  most  aircraft. 
Therefore,  any  signal  processing  technique  that  can  identify  correctly  a  mode  from  2  cycles  of  its 
response  is  a  candidate  technique  for  identifying  the  parameters  of  an  accelerated  aeroelastic  system  in 
a  window  by  window  approach.  Furthermore,  the  analytical  study  typical  wing  section  shows  that 
realistic  level  flight  accelerations  do  not  affect  the  frequencies  and  damping  ratios  of  an  aero  elastic 
system.  In  other  words,  the  aeroelastic  parameters  of  an  aircraft  identified  during  an  accelerated  flight 
are  the  same  as  those  identified  in  stabilized  flight  conditions.  Therefore  flutter  testing  could  be 
performed  in  accelerated  flight,  thereby  reducing  the  cost  and  risk  involved  in  determining  the  flutter 
envelopes  of  fighters.  In  practice,  technical  details  such  as  re-excitation  procedures,  real-time 
identification,  and  flutter  early  alert  systems  must  be  addressed  to  enable  such  a  technology. 
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Therefore  it  can  be  concluded  that  windowing  can  be  used  effectively  to  identify  the 
parameters  using  data  from  an  accelerating  aircraft. 

The  third  issue  was  addressed  in  chapter  7,  namely,  the  need  to  anticipate  the  vibration  response 
at  airspeeds  (or  Mach  numbers)  beyond  which  the  aircraft  has  accelerated.  One  key  to  safely  using  the 
accelerated  flutter  test  method  is  the  development  of  a  predictive  tool  that  can  determine  stability  or  the 
flutter  speed  of  the  system  using  subcritical  flight  test  data.  This  chapter  introduced  the  one  degree  of 
freedom  system  and  explored  different  approaches  to  developing  a  predictive  tool.  The  first  approach 
was  to  develop  an  estimator  that  would  conservatively  predict  the  angle  of  attack  of  the  one  degree  of 
freedom  system.  This  estimator  was  found  to  be  overly  conservative,  and  of  no  practical  use.  The  next 
approach  was  to  find  a  direct  solution  of  the  one  degree  of  freedom  system.  This  approach  naturally 
produced  a  perfect  fit  of  the  one  degree  of  freedom  system,  and  also  did  an  excellent  job  of  predicting 
the  result  of  a  simulation  of  the  two  degree  of  freedom  typical  wing  section.  However  when 
considering  the  stability  of  a  three  dimensional  system  the  simple  predictor  proves  to  be  less  than 
useful.  This  is  an  area  where  further  research  is  needed.  Among  the  potential  solutions  to  this  problem 
are  three  dimensional  analytical  models,  reduced  order  models,  or  full  three  dimensional  simulations 
could  be  used  as  the  predictive  tool.  Ideally  whatever  predictive  tool  is  used  it  should  be  updated  using 
flight  test  data  as  it  becomes  available.  Although  this  issue  merits  more  research  before  it  is  used  for 
live  flutter  testing,  it  has  been  demonstrated  for  a  simple  system  that  a  simplified  model  can  be  used  to 
predict  flutter  parameters  conservatively. 

Therefore  it  can  be  concluded,  at  least  in  principal,  a  simplified  model  of  an  aeroelastic 
system  can  be  used  to  conservatively  anticipate  the  damping  at  airspeeds  higher  than  the 
currently  attained  airspeed. 

The  final  issue  was  also  addressed  in  chapter  7,  that  is  developing  formal  criteria  to  determine 
when  the  aircraft  needs  to  end  the  acceleration  and  terminate  the  test  point.  These  criteria  must  allow 
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some  time  for  aircraft  excitation  (if  decay  data  is  to  be  used)  data  to  be  collected  and  analyzed  and 
corrective  action  taken.  The  biggest  contributors  to  the  time  required  for  this  chain  of  events  are 
aircraft  excitation  which  takes  about  1  second,  data  collection  which  could  take  up  to  one  second  for 
noisy  data,  and  reaction  times  for  the  engineer  and  the  pilot  which  takes  about  one  half  second  each. 
Although  this  three  seconds  could  be  reduced  to  one  second  if  concurrent  excitation  and  analysis  is 
used  and  the  engineer  and  pilot  are  taken  out  of  the  loop,  it  was  shown  in  this  chapter  that  using  the 
three  second  criterion  the  aircraft  could  still  be  halted  in  time  to  avert  a  flutter  situation.  The  Knock-It- 
Off  criteria  must  also  account  for  the  uncertainty  in  the  flight  test  data,  and  analysis  methods,  as  well 
as  the  uncertainty  results  from  the  accelerating  flight  test  procedure.  This  chapter  also  proposed  a 
method  of  handling  the  uncertainties  associated  with  open  air  flight  testing  in  a  conservative  way. 

Therefore  it  can  be  concluded,  given  that  damping  can  be  conservatively  anticipated  for 
higher  airspeeds,  the  aircraft  acceleration  can  be  halted  in  time  to  avoid  entering  a  flutter 
condition. 

After  exploring  each  of  the  issues  above,  it  still  appears  that  the  accelerated  flutter  test  technique 
is  a  viable  approach  to  conducting  flutter  envelope  expansion  flight  testing.  There  are  still  several 
research  issues  that  need  to  be  explored  before  this  method  is  used  to  support  initial  flutter  testing  for  a 
new  aircraft.  In  the  next  chapter  some  of  the  supporting  research  that  need  to  be  accomplished  are 
presented. 
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9.0  Future  Work 


9.1  Improved  Damping  Predictor/Estimator 

One  of  the  keys  to  successfully  implementing  the  accelerated  flutter  testing  capability  is  to  have  a 
reliable  predictor  for  the  damping  values.  For  the  first  implementation  the  results  from  the  full  three- 
dimensional  model  could  be  used.  However  this  approach  provides  predictive  data  for  the  test  points 
as  they  are  planned,  but  not  necessarily  how  they  are  flown.  These  slight  differences  in  altitude,  or 
center  of  gravity  location  add  to  the  uncertainty  that  already  exists.  The  ideal  predictor/estimator 
would  be  capable  of  providing  predicted  damping  values  for  the  aircraft  in  the  configuration  and 
conditions  as  flown  in  near  real  time.  One  approach  to  accomplishing  this  is  through  the  use  of 
Reduced  Order  Modeling  which  is  discussed  in  section  9.2. 


9.2  Reduced  Order  Modeling 

Ideally  the  high  fidelity  3  dimensional  model  would  be  used  to  run  real  time  simulations  to 
support  flight  test  decisions  with  near  real  time  estimates  for  the  maneuvers  as  flown  as  opposed  to  as 
planned.  However  existing  compute  resources  are  not  adequate  to  run  the  model  anywhere  near  real 
time.  It  is  still  desirable  to  have  reliable  results  to  compare  to  flight  test  data.  One  method  of  reaching 
these  results  in  a  timely  manner  is  to  use  a  reduced  order  model.  In  order  to  develop  the  reduced  order 
model  the  results  of  the  full  order  model  are  needed.  The  algorithm  uses  these  results  to  estimate  a 
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reduced  order  system  that  produces  the  same  observable  dynamics  within  the  allowable  error  bound  as 
the  full  order  model.  The  reduced  order  model  is  computationally  much  less  expensive  and  will 
therefore  run  very  near  real  time. 


9.3  Flight  Test  Data  Based  Predictor/Estimator 

Even  though  the  full  order  model  is  very  high  fidelity,  and  therefore  can  do  an  excellent  job  of 
predicting  the  dynamics  of  the  aircraft  as  designed,  or  as  represented  in  the  model,  there  will  always  be 
discrepancies  between  the  model  as  designed  (or  modeled)  and  the  real  article  as  built.  Because  of  this 
discrepancy,  the  simulation  data  is  always  somewhat  suspect  until  it  is  validated  with  flight  test  data. 
However,  occasionally  the  validation  process  will  only  reveal  the  discrepancy  between  the  model  and 
the  real  article.  It  is  in  this  case  that  it  is  desirable  to  update  the  model  such  that  it  more  accurately 
reflects  the  dynamics  observed  through  flight  test.  Regardless  of  whether  the  model  is  the  full  order 
model  or  the  reduced  order  model,  if  it  is  being  used  to  predict  or  estimate  the  behavior  of  the  aircraft 
at  higher  Mach  numbers  than  has  been  achieved  through  flight  test,  it  is  critical  that  the  model  reflect 
the  real  life  dynamics. 


9.4  Evaluate  Other  Parameter  Identification  Methods 

Eigensystem  Realization  Algorithm  was  used  exclusively  to  perform  parameter  identification  in 
this  thesis.  There  are  many  other  state  space  parameter  identification  algorithms  that  may  perform  as 
well  or  better  than  ERA  for  the  aircraft  flutter  application.  Although  ERA  has  met  the  initial 
requirements  it  would  be  prudent  to  extend  the  investigation  to  other  parameter  identification  methods. 
Areas  in  which  ERA  could  be  significantly  improved  upon  are, 


1) 


Concurrent  excitation  and  analysis 
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2)  Forgetting  factor 

3)  Noise  tolerance 

Concurrent  excitation  and  analysis  enable  the  flight  test  engineer  to  use  data  that  is  collected 
during  the  excitation  of  the  aircraft  to  perform  the  parameter  identification.  This  would  shorten  the 
time  required  to  look  ahead  and  would  also  provide  nearly  continuous  damping  estimations  as  the 
aircraft  accelerates.  Both  of  these  effects  would  improve  the  safety  of  conducting  accelerating  flutter 
testing. 

A  forgetting  factor  weights  the  data  points  at  the  end  of  the  window  more  heavily  than  the  data 
points  at  the  beginning  of  the  window.  The  weighting  should  allow  the  parameters  to  be  identified 
more  closely  to  the  end  of  the  window  which  would  reduce  the  lag  time  between  the  time  the 
parameter  represents  and  the  time  the  result  is  available  to  make  decisions. 

ERA  performed  reasonably  with  regard  to  noise  tolerance,  however  it  required  five  times  as  much 
data  (0.2  seconds  for  clean  data  vs.  1 .0  seconds  for  noisy  data)  to  identify  parameters  from  noisy  data. 
The  ability  to  identify  parameters  using  a  smaller  window  of  data  will  improve  the  safety,  and/or  the 
efficiency  of  the  accelerated  flutter  test  technique. 


9,5  Refine  Knoek-It-Off  Criteria 

The  Knock-It-Off  criteria  presented  in  chapter  7  contains  the  essential  characteristics  that  are 
necessary  to  support  initial  flutter  testing  on  a  new  aircraft.  These  criteria  provide  a  good  outline  but 
need  to  be  further  developed  by  soliciting  requirements  from  the  flutter  flight  test  community.  Clearly 
the  final  Knock-It-Off  criteria  should  incorporate  the  improved  damping  predictor/estimator  discussed 
above.  It  would  also  be  helpful  to  develop  statistical  basis  for  evaluating  the  effect  that  noise  has  on 
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the  uncertainty  of  the  damping  value  identified  from  the  flight  test  data.  However,  it  may  be  difficult 
to  attain  enough  data  to  do  any  kind  of  statistical  analysis. 
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